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ABSTRACT 


Two  numerical  methods  for  solving  the  incompressible  Navier-Stokes  equations  are 
compared  with  each  other  by  applying  them  to  calculate  laminar  and  turbulent  flows 
through  curved  ducts  of  regular  cross-section.  Detailed  comparisons,  between  the 
computed  solutions  and  experimental  data,  are  carried  out  in  order  to  validate  the  two 
methods  and  to  identify  their  relative  merits  and  disadvantages.  Based  on  the  conclusions 
of  this  comparative  study  a  numerical  method  is  developed  for  simulating  vi.scous  flows 
through  curved  ducts  of  varying  cross-sections.  The  proposed  method  is  capable  of 
simulating  the  near-wall  turbulence  using  fine  computational  meshes  across  the  sublayer  in 

conjunction  with  a  two-layer  k-e  model.  Numerical  solutions  are  obtained  for:  i)  a  straight 
transition  duct  geometry,  and  ii)  a  hydroturbine  draft-tube  configuration  at  model  scale 
Reynolds  number  for  various  inlet  swirl  intensities.  The  report  also  provides  a  detailed 
literature  survey  that  summarizes  all  the  expierimental  and  computational  work  in  the  area  of 
duct  flows. 


ACKNOWLEDGEMENTS 

The  research  described  herein  was  sponsored  by  a  grant  from  the  Tennessee  Valley 
Authority,  with  supplementary  funds  from  the  U.S.  Bureau  of  Reclamation.  The  authors 
are  most  grateful  to  Bill  Waldrop,  Paul  Hopping,  and  Patrick  March,  of  TVA,  for  their 
advice  and  support.  Thanks  are  also  due  to  Tom  McKay,  of  EG&G-ldaho  Inc.,  for  his 
help  in  obtaining  access  to  the  supercomputers  of  the  Idaho  National  Engineering 
Laboratory. 

The  calculations  presented  herein  were  carried  out  on  the  Cray  2  supercomputer  and 
the  IBM  RS-600()  superscalar  computers  of  the  National  Center  for  Supercomputing 
Applications  at  the  University  of  Illinois  at  Urbana-Champaign,  the  CRAY  YMP 
supercomputers  of  the  Numerical  Aerodynamic  Simulation  facility  at  NASA  Ames 
Research  Center  and  the  CRAY  X-MP/216  supercomputer  at  the  Idaho  National 
Engineering  Laboratory. 


Ill 


LIST  OF  SYMBOLS 


Alphabetical  Symbols 


A,  B,  C,  R, 
a,b,c 

A(j),  B(|), 

Ag 

Cd,  Cp,  Cu, 

c, 

Cnb 


coefficients  in  the  linearized  transport  equations  for  the  finite 
analytic  method 

constants  in  the  finite  analytic  solution  of  the  one-dimensional 
equations 

diagonal  matrices  containing  the  convective  velocities  (j  =  1,2,3) 

coefficients  in  the  linearized  transport  equations  for  (|) 

constants  in  the  two-layer  turbulence  model 

finite-analytic  coefficients 

constant  in  two-layer  turbulence  model 

finite-analytic  coefficients 


Cp 


pressure-coefficient  (  =  2(P  -  Po)/  pUo“  ) 


Cgi,  Ce2 

CFL 

D 

De 

dh 

E‘‘ 


f> 

G 

g 


Gp 

g'j 

gij 


turbulence  model  constants 
Courant-Friedrich-Lewis  number 

mass-source  in  the  pressure  equation  of  Method  1 

Dean's  number 

duct's  hydraulic  diameter 

coefficients  in  the  pressure  equation  of  Method  1 

viscous  flux  vector  (i  =  1,2,3) 

grid  control  functions  (i  =  1 ,2,3) 
turbulence  generation  term 

determinant  of  the  metric  tensor  gjj  and  source  function  at  the  node 
P  for  the  finite  analytic  discretization 

source  term  at  the  node  P  for  the  finite-analytic  discretization 

contravariant  metric  tensor 

covariant  metric  tensor 


1  v 


h,  k,  1 

J 

k 

it 

P 

Q 

r 

Re 

R0 

S 

S(|) 

t 

U,  V.  w. 


lib 


Ui 

Uo 


Ut 

UiUj 

V» 

V(i) 

V(i) 

X,  Y,  Z 


xi 


dimensions  of  a  finite-analytic  numerical  element 
Jacobian  of  the  geometric  transformation 
turbulent  kinetic  energy 

length  scales  in  the  two-layer  turbulence  models 

piezometric  pressure,  normalized  by  pUo^ 

velocity  vector  containing  the  three  cartesian  velocity  components 

radial  coordinate 
Reynolds  number 

effective  Reynolds  number  including  the  eddy  viscosity 
inlet  swirl  intensity 

source  functions  in  the  transport  equations,  used  in  Method  1 
dimensionless  time,  normalized  by  Lo/Uo 

dimensionless  mean  velocity  components  in  the  streamwise,  radial 
and  normal  directions  (used  for  curved  duct  geometries)  and 
cartesian  velocity  components  in  the  x-,  y-,  and  z-directions 

bulk  velocity 

dimensionless  mean  Cartesian  velocity  component 
free-stream  (reference)  velocity 
friction  velocity  (=  Vx^/p  ) 
dimensionless  Reynolds  stresses  (i,  j  =  1 ,2,3) 
contravariant  velocity  components  (i  =  1 ,2,3) 
physical  contravariant  velocity  components  (i  =  1 ,2,3) 

pseudo-velocity  of  V(i) 

curvilinear  coordinates  along  the  curved  duct 

cartesian  coordinates  (i  =  1, 2,  3) 

Uxy  ,, 

dimensionless  wall  distance  (  =  ) 


V 


z 

Greek  Symbols 
a( 

Y 

rjk 

5 

At.  At 
e 
0 
K 

V 

Vt 

Ti,i; 

tJ(j) 

X 

Xw 

1) 

tOi 

a 


u^y 

dimensionless  wall  distance  ( =  ) 

normal  coordinate 

Runge-Kutta  coefficients 
artificial  mass  source  parameter 
Cristoffel  symbols  of  the  second  kind 
finite  difference  operator 
grid  spacings  (i  =  1,2,3) 
time  increments 

rate  of  turbulent  energy  dissipation,  normalized  by  U  /  L 

streamwise  coordinate 

Von  Karman  constant  (  =  0.42  ) 

kinematic  viscosity 

turbulent  eddy-viscosity 

general  curvilinear  coordinates 

general  curvilinear  coordinates  ( i  =  1,  2.  3  ) 

turbulence  model  constants  ((])  =  U,  V,  W,  k,  e) 

time  variable 

wall  shear  stress 

transport  quantities  and  circumferential  coordinate 
implicit  residual  smoothing  coefficients  (i  =  1,2.3) 

Von- Neumann  stability  number 
streamwise  mean  vorticity  component 


vi 


I.  INTRODUCTION 


The  study  of  flow  through  curved  ducts,  or  conduits  is  of  great  importance  because  such 
flows  are  found  in  a  variety  of  practical  situations,  ranging  in  scale  from  the  flow  through  the 
human  arterial  and  respiratory  systems  to  that  in  the  ducting  of  gas  and  hydraulic  turbines.  The 
flow  in  open  channels  and  rivers  is  also  related  to  duct  flow,  although,  in  these  cases,  both  the 
bottom  and  the  free  surface  may  be  free  to  move.  An  accurate  description  of  the  flow  in  ducts, 
even  when  the  walls  are  rigid,  poses  a  challenging  fluid  mechanics  problem.  There  are  several 
difficulties.  One  of  these  is  the  complexity  of  the  geometries  encountered  in  practice.  Changes  in 
cross-sectional  shape  and  curvature  of  the  duct  axis  are  of  particular  concern  because  these  lead  to 
vortical  motions,  flow  reversals,  and  unsteadiness.  Another  major  difficulty  is  related  to  the  fact 
that  most  flows  of  practical  interest  are  turbulent,  and  available  mathematical  models  cannot 
properly  simulate  all  of  the  consequences  of  turbulence.  To  the.se  may  be  added  difficul-ies  that  are 
peculiar  to  specific  applications,  such  as  coupling  between  the  flow  and  moving  boundaries  in 
biomedical  applications,  compressibility  effects  and  shock  waves  in  gas  turbines,  sediment 
transport  in  rivers,  and  aeration  in  autoventing  hydraulic  turbines.  These,  and  other  applications 
too  numerous  to  list,  have  inspired  many  experimental,  analytical  and  computational  studies  of  the 
flow  in  curved  ducts  but,  for  a  variety  of  reasons,  it  is  .still  not  possible  to  accurately  predict  such 
flows. 

The  work  described  in  this  report  was  motivated  by  a  recently  initiated  joint  research 
program  between  the  Iowa  Institute  of  Hydraulic  Research  (IIHR)  and  the  Tennes.see  Valley 
Authority  (TVA).  This  program  aims  to  develop  a  numerical  method  specially  for  applications  to 
the  new  generation  of  autoventing  hydroturbines  (AVT)  that  are  being  considered  as  possible 
replacements  for  conventional  ones.  The  initial  phase  of  the  project  is  conducted  in  two  parts: 
development  of  a  numerical  method  for  solution  of  the  Reynold.s-averaged  Navier-Stokes  (RANS) 
equations  for  single-phase  flow  in  ducts  of  complex  geometry  is  pursued  at  IIHR.  and 
development  of  physical  models  and  transport  equations  for  two-phase  flow,  for  air-water 
exchange,  is  carried  out  at  TVA.  The  two  phases  will  be  integrated  at  an  appropriate  stage.  This 
report  describes  the  development  of  numerical  methods  for  prediction  of  single-phase  fluid  flow  in 
curved  ducts. 

The  starting  point  for  the  research  described  here  was  provided  by  the  two  numerical 
methods  that  were  available  at  IIHR  at  the  beginning  of  the  project,  namely,  one  developed  by 
Patel  and  his  colleagues  (.see,  for  example,  Chen,  Patel  and  Ju  (1990)),  and  the  other  developed  by 
Sotiropoulos  (1991).  Both  were  developed  for  external,  three-dimensional,  laminar  and  turbulent 
flows,  and  have  been  extensively  applied  to  such  flows.  Both  methods  were  modified  to 
incorporate  the  two-layer  turbulence  model  of  Chen  and  Patel  (1988),  and  extended  to  internal  flow 


in  ducts.  As  a  preliminary  step  in  the  validation  of  these  methods,  consideration  is  first  given  to 
the  flow  in  simply-curved  ducts  of  simple  cross-section.  The  choice  of  suitable  test  cases  to 
evaluate  the  performance  of  these  methods  is  not  straightforward,  however,  because  there  exists  an 
enormous  amount  of  previous  work  on  this  subject  and  even  a  cursory  review  of  this  indicates  that 
there  are  a  number  of  issues,  both  numerical  and  physical,  that  are  not  yet  fully  resolved. 
Therefore,  this  report  begins  with  a  brief  review  of  previous  studies  of  curved  duct  flows  with 
emphasis  on  experiments  and  those  aspects  that  are  still  topics  of  active  research.  The  two 
methods  are  then  described,  and  their  results  are  compared  for  a  few  relatively  simple  duct  flows  to 
highlight  their  capabilities  and  shortcomings.  Next,  a  modified  version  of  the  method  of 
Sotiropoulos  is  employed  to  study  the  flow  in  additional  cases.  Among  these  are  the  flow  in  a 
straight  duct  that  changes  its  cross  section  from  circular  to  rectangular,  and  the  flow  in  a  typical 
draft  tube  of  a  hydraulic  turbine.  These  two  examples  bring  forth  the  complexities  that  are  present 
when  the  duct  cross-section  shape  is  changed  and  the  area  is  increased  for  continuous  diffusion. 
The  results  demonstrate  that  accurate  numerical  simulation  of  such  practical  flows  is  at  hand 
provided  some  improvements  are  made  in  two  areas,  namely,  in  the  turbulence  model  and  in 
generation  of  grids  appropriate  for  more  complex  geometries.  Possible  future  directions  of  the 
project  are  outlined  in  the  final  section  of  the  report. 

11.  OVERVIEW  OF  PREVIOUS  WORK  ON  FLOW  IN  CURVED  DUCTS 

Before  attempting  to  summarize  previous  work  on  flow  in  curved  ducts,  it  is  useful  to 
define  certain  terms  that  are  commonly  used  in  connection  with  such  flows.  Definition  of  these 
terms  aids  in  the  classification  of  duct  flows,  in  general,  and  provides  a  basis  for  comparing  the 
various  types  of  flows  that  have  been  studied. 

The  flow  in  a  curved  duct  of  finite  cross-section  is  three  dimensional  insofar  as  the  velocity 
vector  has  three  non-zero  components,  each  of  which  varies  across  the  section.  Whether  the  flow 
is  laminar  or  turbulent,  such  a  three-dimensional  flow  is  often  described  in  terms  of  primary  and 
secondary  velocities.  Loosely  speaking,  primary  refers  to  the  velocity  component  along  the  axis  of 
the  duct  and  secondary  flow  refers  to  velocity  components  in  transverse  sections,  orthogonal  to  the 
duct  axis.  For  ducts  whose  cross-section  shape  and  area  change  rapidly,  these  terms  are  not 
precise  because  these  velocity  components  are  not  orthogonal  and,  therefore,  depend  on  the 
coordinates  that  are  chosen.  The  term  secondary  also  suggests  that  motion  is  weak  compared  with 
the  primary  motion.  This  is  again  not  true  for  ducts  of  arbitrary  shape  and  curvature.  However, 
most  previous  studies  of  curved-duct  flow  are  confined  to  simple  cross  sections  and  moderate 
curvatures  for  which  the  concepts  of  primary  and  secondary  flow  are  quite  useful. 
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In  a  curved  duct,  the  secondary  flow  arises  quite  simply  as  a  result  of  the  radial  pressure 
gradient.  The  higher  pressure  at  the  outer  radii  drives  the  slower  moving  near-wall  fluid  towards 
the  inner  radii.  This  is  compensated  by  ai:  opposite  but  considerably  weaker  motion  of  the  higher- 
velocity  fluid  in  the  duct  core  from  the  inner  to  the  outer  radii.  Prandtl  termed  the  resulting 
secondary  flow  as  secondary  motion  of  the  first  kind  or  pressure-driven  secondary  motion.  In 
turbulent  flow,  another  kind  of  secondary  motion  can  be  induced,  even  in  straight  ducts,  by  the 
gradients  of  the  Reynolds  stresses  in  the  plane  normal  to  the  primary  flow.  Prandtl  called  this 
secondary  motion  of  the  second  kind.  This  stress-driven  secondary  motion  is  absent  in  laminar 
flow,  but  in  turbulent  flow,  both  are  present.  For  an  overview  of  the  theory  of  secondary  motion 
in  turbulent  flows,  reference  should  be  made  to  the  recent  article  by  Bradshaw  (1987). 

The  magnitude  of  the  pressure-driven  .secondary  motion  depends  on  a  number  of  factors, 
including  the  section  shape,  the  Reynolds  number,  the  Dean  number,  and  the  flow  conditions  at  the 
beginning  of  curvature.  In  general,  Reynolds-stress  driven  secondary  motion  is  much  weaker  than 
pressure-driven  secondary  motion  but  the  two  are  not  easily  separated  because  of  the  coupling 
between  the  mean  flow  and  the  Reynolds  stresses.  Secondary  motions  of  both  kind  have  a 
significant  impact  on  the  distribution  of  the  wall  .shear,  heat  transfer  at  the  wall,  and  related 
properties. 

A  review  of  previous  work  on  curved  ducts  leads  to  the  following  general  observations. 
Most  studies  have  dealt  with  ducts  of  simple  shape,  such  as  a  circular  or  square  section,  with 
constant  cross-sectional  area,  turning  through  either  90  or  180  degrees  with  a  constant  radius. 
Ducts  with  two  bends,  the  so-called  S-shaped  ducts,  also  fall  in  this  category.  There  are  very  few 
investigations  in  ducts  with  section  and  area  changes,  or  in  ducts  with  axis  curved  in  more  than  one 
plane.  Secondly,  there  are  many  more  studies  of  laminar  flow,  at  necessarily  low  Reynolds 
numbers,  than  with  turbulent  flow.  Transitional  flows  have  not  been  investigated.  Thirdly,  fully- 
developed  entry  flow  is  more  common  than  developing  boundary-layer  flow  at  entry.  The  former 
leads  to  rapidly  developing  and  stronger  secondary  motions.  Fourth,  in  spite  of  extensive 
experimental  investigations,  the  data  base  is  still  limited,  particularly  with  regard  to  the  distribution 
of  the  Reynolds  stresses  in  turbulent  flow. 

The  better  documented  experiments  in  laminar  and  turbulent  flows  have  been  used  to  test 
computations  using  a  variety  of  numerical  methods  and  turbulence  models.  Of  course,  there  are 
many  computations  for  each  experimental  condition,  and  the  results  have  been  rather  mixed. 
Differences  among  numerical  methods  and  turbulence  models  will  be  di.scussed  in  the  text  as  the 
results  of  the  present  methods  are  presented  and  compared  with  previous  ones.  For  the  present 
purposes,  it  is  more  useful  to  provide  a  summary  of  some  of  the  most  commonly  cited 
experiments.  This  summary  is  not  intended  to  be  exhaustive,  however.  We  have  .selected  the  most 
recent  experiments,  performed  with  advanced  instrumentation  to  obtain  the  most  detailed  data. 
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II.  1  Experiments  in  laminar  flow 

Laminar  flow  through  curved  pipes  have  been  studied  extensively  during  the  past  fifteen 
years,  primarily  because  of  their  relevance  to  biomedical  problems.  Measurements  for  steady  and 
unsteady  flows  have  been  reported  by,  among  others,  Agrawal  et  al.  (1978),  Chandran  et  al. 
(1979,  1981),  Enayet  et  al.  (1982b),  Talbot  and  Gong  (1983),  Bovendeerd  et  al.  (1987)  and  Rindt 
et  al.  (1991).  There  are  several  other  experiments  in  curved  ducts  of  different  shapes  but  none 
provide  data  thci  are  comparable  in  detail  to  the  experiments  summarized  in  Table  1.  Together, 
these  experiments  provide  a  data  base  for  the  development  of  secondary  motion  in  laminar  flow 
with  developed  and  developing  flow  at  entry. 

II. 2  Reynolds-stress  driven  secondary  motion:  experiments  in  turbulent  flow 
in  straight  ducts 

Secondary  motions  of  the  second  kind  were  first  observed  by  Nikuradse  ( 1926)  who  made 
careful  measurements  of  mean  velocity  in  turbulent  flow  in  a  series  of  straight,  non-circular  ducts. 
His  measurements  indicated  that  the  isovels  (or  isotachs,  contours  of  constant  axial  velocity)  were 
distorted  near  the  comer  region  in  a  manner  that  suggested  the  existence  of  secondary  motion. 
This  situation  is  depicted  in  Figure  1,  taken  from  Gessner  and  Jones  (1965),  where  a  typical  isovel 
pattern  for  laminar  flow  is  compared  with  the  corresponding  one  for  turbulent  flow.  Recall  that 
there  is  no  secondary  motion  in  laminar  flow.  This  has  been  shown  analytic:'llv  by  Moissis  ( 1957) 
and  Maslen  (1958).  It  is  seen  that  in  the  turbulent  case  the  isovels  are  displaced  towards  the  comer 
and  away  from  the  mid-point  of  the  walls.  Prandtl  (1926)  suggested  that  the  distortion  of  the 
isovels  is  the  result  of  a  secondary  motion  towards  the  corner,  which  is  accompanied  by  a  return 
flow  at  the  mid-points  of  the  walls  for  continuity  to  be  satisfied.  Moreover,  he  postulated  that  this 
secondary  motion  is  generated  by  velocity  fluctuations  tangential  to  the  isovels  in  regions  where  the 
curvature  of  the  isovel  changes. 

Due  to  the  small  magnitude  of  the  secondary-flow  velocity  components,  approximately  one 
percent  of  the  average  or  bulk  velocity,  quantitative  measurements  for  developing  turbulent  flows 
in  square  ducts  were  not  reported  until  Hoagland  (1960)  developed  a  hot-wire  technique. 
Subsequently,  hot-wire  measurements  of  increased  accuracy  were  reported  by  Leutheusser  (1963), 
Brundrett  and  Baines  (1964),  Gessner  (1964),  Gessner  and  Jones  (1965),  Ahmed  (1971),  Thomas 
and  Easter  (1972),  and  Launder  and  Ying  (1972).  The  first  attempt  to  provide  a  comprehensive 
explanation  of  the  physics  of  stress-driven  secondary  flow  was  that  of  Brundrett  and  Baines  who 
measured  the  three  mean-velocity  components  and  the  six  Reynolds-stress  components.  In  their 
analysis,  they  employed  the  equation  for  the  mean  streamwise  vorticity  along  with  the  symmetry 
properties  the  Reynolds-stress  tensor,  and  showed  that:  (i)  streamwise  vorticity  is  produced  in  the 
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region  near  the  corner  bisector  by  the  gradients  of  the  Reynolds  stresses,  (ii)  the  corner  bisectors 
separate  independent  secondary  flow  circulation  zones  (as  seen  in  figure  lb),  and  (iii)  the  main 
contribution  to  streamwise  vorticity  production  comes  from  the  gradients  of  the  normal  Reynolds 
stresses.  The  validity  of  their  last  conclusion,  however,  was  later  challenged  by  Perkins  (1970) 
and  Gessner  (1973).  Perkins,  for  instance,  pointed  out  that  the  experimental  apparatus  of  Brundett 
and  Baines  could  lead  to  ±100%  error  in  the  shear  stress  measurement,  while  Gessner  (1973) 
claimed  that  the  gradients  of  the  normal  Reynolds  stresses  do  not  play  a  major  role  in  the 
streamwise  vorticity  production  and  it  is  rather  the  transverse  gradients  of  the  shear  stress  that  drive 
the  secondary  motion  in  the  corner  region.  In  a  more  recent  work,  Demuren  and  Rodi  (l'^84) 
reviewed  all  the  above  experimental  findings  and  concluded  that  both  the  gradients  of  normal  and 
shear  Reynolds  stresses  are  of  the  same  order  of  magnitude  and  that  it  is  their  difference  that  drives 
the  secondary  motion. 

In  most  of  the  above  experimental  efforts  ihe  emphasis  was  placed  in  understanding  the 
physical  mechanism  that  generates  and  sustains  the  secondary  motion,  rather  than  in  providing  a 
reliable  data  base  suitable  for  validating  turbulence  models.  In  many  of  these  studies,  for  instance, 
the  symmetry  of  the  flow  about  the  corner  bisector  was  either  not  checked  (since  measurements 
were  carried  out  in  only  one  octant),  or  the  check  was  not  thorough  enough  to  guarantee  reliability 
of  the  measurements.  In  a  more  recent  study,  Melling  and  Whitelaw  (1976)  attempted  to  obtain 
better  quality  measurements  in  the  entrance  region  of  square  ducts  using  Laser-Doppler 
Velocimetry  (LDV).  They  made  measurements  in  a  quadrant  of  the  flow  at  several  cross-sections 
and  reported  high  levels  of  symmetry  in  the  streamwise  mean  velocity  and  normal  Reynolds-stress 
components  (the  maximum  asymmetry  was  found  to  be  about  3  percent  in  the  near  wall  region  and 
1  percent  elsewhere).  The  same  degree  of  symmetry,  however,  was  not  achieved  in  their 
secondary  flow  velocities  and  transverse  Reynolds-stress  components,  and  their  data  indicate 
significant  asymmetries  about  the  wall  bisector.  Another  set  of  measurements  in  developing 
turbulent  flow  in  the  entrance  region  of  a  square  duct  was  reported  by  Gessner  et  al.  ( 1977);  see 
also  Po  (1975).  They  carried  out  detailed  mean  velocity  (in  a  quadrant  of  the  flow)  and  turbulence 
(in  an  octant  of  the  flow)  measurements  at  several  cross-sections  in  the  developing  and  fully- 
developed  flow  regions.  Their  mean  velocity  data  exhibit  high  levels  of  symmetry,  but  the  quality 
of  their  turbulence  data  can  not  be  readily  evaluated.  The  measurements  of  Gessner  et  al.  constitute 
a  complete  and  comprehensive  data  set  which  is  detailed  enough  to  assist  in  the  development  and 
validation  of  turbulence  models.  A  detailed  review  of  experiments  and  computations  for  straight 
ducts  of  regular  cros.s-section  can  be  found  in  Demuren  and  Rodi  (1984).  Additional  experimental 
data  for  duct.'’  with  arbitrary  cross-sections  have  been  reported  by  Rodet  (1960),  Aly  et  al.  (1978) 
and  Seale  (1982)  for  trapezoidal,  triangular,  and  simulated  rod-bundle-type  cross-sections, 
respectively. 
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As  already  noted,  the  driving  force  for  secondary  motions  of  the  second  kind  is  the 
anisotropy  of  the  Reynolds  stresses  and,  thus,  any  attempt  to  calculate  such  motions  should  utilize 
a  turbulence  model  which  can  predict  turbulence  anisotropy.  It  has  been  demonstrated 
computationally  (see  Demuren  and  Rodi,  1984)  that  the  isotropic  k-e  model  does  not  predict  any 
secondary  motion  of  second  kind  and  fails,  consequently,  to  predict  the  important  features  of 
straight  duct  flows.  For  that  reason  all  the  computational  efforts  in  the  area  of  stress-driven 
secondary  motions  have  focused  on  developing  algebraic  Reynolds  stress  (ASM),  or  full  Reynolds 
stress  (RSM)  models.  The  first  calculation  of  secondary  motion  of  the  second  kind,  in  straight 
non-circular  ducts,  was  reported  by  Launder  and  Ying  (1973).  They  employed  an  ASM  model 
(LY  model,  ^ay)  based  on  a  simplification  of  the  model  of  HanjaJic  and  Launder  (HL)  (1972).  The 
LY  model,  or  some  variant  of  it,  was  subsequendy  used  by  Gessner  and  Emery  (1981),  Trupp  and 
Aly  (1979),  Rapley  and  Gosman  (1986),  and  Nakayama  et  al.  (1983)  to  calculate  flows  through 
straight  ducts  of  regular  and  irregular  cross-sections.  All  the  above  calculations  yielded  fairly  good 
results  for  the  magnitude  of  the  secondary  flow  but  failed  to  yield  accurate  predictions  for  the 
Reynolds  stress  distributions.  Naot  and  Rodi  (NR)  (1982)  were  able  to  obtain  better  overall 
predictions,  for  both  the  crossflow  and  the  Reynolds  stress  distributions,  by  utilizing  an  ASM 
model  which  they  developed  from  the  Reynolds  stress  model  of  Launder,  Reece  and  Rodi  (LRR) 
(1975).  A  more  refined  version  of  the  NR  model  was  proposed  by  Demuren  and  Rodi  (DR) 
(1984,  1987)  and  used  to  calculate  the  flow  in  ducts  with  simple  cross-sections.  Recently, 
Demuren  (1991)  has  extended  both  the  NR  and  DR  models  to  generalized  curvilinear  coordinates 
and  applied  them  to  calculate  flows  in  straight  ducts  of  arbitrary  cross-sections.  Most  of  his 
computations  demonstrate  correct  trends  but  are,  at  best,  in  qualitative  agreement  with  the 
experimental  data.  It  should  be  pointed  out  that  all  of  the  above  described  computational  methods 
utilize  the  wall-function  approach  which,  as  indicated  by  the  limited  success  of  even  the  most 
advanced  methods,  is  inadequate  for  complex  shear  flows  (Demuren,  1991). 

II. 3  Experiments  in  turbulent  flow  in  curved  ducts  of  regular  cross- 

sectional  shape 

Turbulent  flow  measurements  in  curved  ducts  of  square  or  rectangular  cross-section  have 
been  reported  by  Joy  (1950),  Eichenberger  (1952,  1953),  Squire  (1954),  Bruun  (1979), 
Humphrey  et  al.  (1981),  Taylor  et  al.  (1982),  Enayet  et  al.  (1982a),  Chang  (1983),  lacovides  et  al. 
(1990)  and  Kim  (1991).  The  early  works  of  Joy,  Eichenberger  and  Squire  provided  evidence  of 
the  oscillatory  nature  of  the  secondary  flow  in  180^  square  bends.  The  turbulent  flow 
measurements  of  Joy,  for  instance,  indicate  that  the  secondary  flow  was  reversed  at  a  bend  angle 
of  135°.  Bruun  (1979)  studied  the  effect  of  inlet  boundary-layer  thickness  on  the  development  of 
the  secondary  flow  and  on  the  pressure  losses  in  a  120°  square  bend.  He  presented  data  for  the 
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mean  total  pressure,  static  pressure  and  velocity  field  as  well  as  some  limited  turbulence  intensity 
data. 

Table  2  summarizes  the  most  detailed  experiments  in  turbulent  flow  in  curved  ducts.  A 
significant  contribution  in  the  understanding  of  the  physics  of  curved-duct  flow  came  from 
Humphrey  et  al.  (1981)  who  presented  detailed  measurements  inside  a  highly  curved  90^^  square 
bend.  They  employed  LDV  to  measure  the  mean-velocity  components  and  components  of  the 
Reynolds-stress  tensor  at  several  cross  sections  inside  the  bend.  Their  data  indicates  that  the 
Reynolds-stress  driven  secondary  motion,  measured  in  the  inlet  tangent  of  the  bend,  was  quickly 
taken  over  by  the  pressure-driven  motion,  in  which  the  secondary  velocities  attain  values  up  to  28 
percent  of  the  mean  bulk  velocity,  during  the  first  half  of  the  bend.  This  latter  secondary  motion 
was  found  to  be  responsible  for  strong  cross-stream  convection  of  the  Reynolds  stresses  inside  the 
bend  which,  in  conjunction  with  the  curvature  effects,  resulted  in  a  high  level  of  anisotropy  of  the 
turbulence.  Furthermore,  the  production  of  turbulent  kinetic  energy  predominated  near  the  outer 
wall  (destabilization  due  to  concave  curvature)  but  regions  of  negative  production  were  observed 
near  the  inner  wall  (stabilization  by  convex  curvature).  To  study  the  effect  of  the  inlet  boundary- 
layer  thickness  on  the  flow  structure,  Taylor  et  al.  (1982)  carried  out  meaiurcments  in  exactly  the 
same  bend  but  with  developing  flow,  compared  to  the  nearly  fully-developed  flow  in  the 
experiments  of  Humphrey  et  al.,  at  the  inlet  to  the  bend.  They  concluded  that  the  thinner  inlet 
boundary  layer  results  in  weaker  crossflow  and  lower  overall  turbulence  intensities. 

The  effect  of  bend  curvature  on  the  crossflow  and  the  turbulence  structure  was  studied  by 
Enayet  et  al.  (1982a)  who  carried  out  measurements  inside  a  90*^  bend  of  milder  curvature  with  a 
partially  developed  flow  (relatively  thick  boundary  layer)  at  inlet.  Their  data  shows  that,  due  to  die 
milder  curvature  the  streamwise  pressure  gradient  is  small  and  less  important.  Also,  unlike 
observations  in  strongly  curved  ducts,  there  appears  to  be  no  evidence  of  significant  damping  or 
amplification  of  turbulence  due  to  longitudinal  curvature.  Taylor  et  al.  compared  their 
measurements  with  the  data  of  Enayet  et  al.  and  concluded  that  the  milder  curvature  results  in  much 
weaker  crossflow  velocities  which  are  confined  to  a  smaller  region  near  the  wall.  However,  due  to 
the  small  streamwise  pressure  gradient,  in  the  case  of  mild  curvature,  the  effect  of  the  crossflow  in 
the  streamwise  flow  development  was  found  to  be  dominant. 

The  flow  through  a  strongly  curved  180”  square  bend  was  studied  by  Chang  (1983).  His 
detailed  measurements  (wall  pressures,  mean  velocity  components  and  Reynolds  stresses)  revealed 
rather  complicated  flow  patterns,  with  the  secondary  motion  dominating  the  dynamics  of  the  flow 
development.  Between  0”  and  45”,  for  instance,  destabilizing  and  stabilizing  curvature  effects  in 
conjunction  with  the  strong  secondary  motion-primarily  pressure  driven  except  in  the  inlet  straight 
tangent,  where  Reynolds  stress  driven  secondary  currents  are  present— cause  large  levels  of 
turbulence  anisotropy.  As  a  result,  a  region  of  negative  production  of  turbulent  kinetic  energy 
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exists  near  the  inner  wall,  indicating  flow  of  energy  from  the  large  eddies  back  to  the  mean  flow 
(recall  that  a  similar  trend  was  observed  by  Humphrey  et  al.  (1981)  in  their  90°  bend  experiment). 
The  anisotropy  of  the  turbulence  field  is  reduced,  however,  for  bend  angles  greater  than  130°,  due 
to  the  mixing  effect  of  the  crossflow.  Apart  from  influencing  the  dynamics  of  turbulence,  the 
secondary  motion  has  a  significant  impact  on  the  mean  flow  development  as  well.  For  bend  angles 
less  than  90°  the  strong  secondary  currents  convect  low  momentum  fluid  from  the  side  walls  to 
the  core  of  the  cross-section,  resulting  in  S-shaped  axial  and  radial  mean  velocity  profiles.  The 
same  phenomenon  was  observed  by  Humphrey  et  al.  ( 1 977)  in  their  laminar  flow  measurements, 
but  it  was  not  present  in  the  turbulent  flow  measurements  of  Humphrey  et  al.  (1981)  despite  the 
fact  that  the  same  bend  was  used  in  both  experiments.  For  angles  greater  than  90°,  however,  the 
crossflow  starts  decreasing  and  high  speed  fluid  is  restored  back  in  the  core  of  the  flow.  In  a 
more  recent  study,  lacovides  et  al.  (1990)  carried  out  turbulent  flow  measurements  through  a  180° 
square  bend  whose  geometry  was  similar  to  that  used  by  Ghang  ( 1983)  except  that  the  inlet  section 
was  shortened  to  six  hydraulic  diameters  (as  compared  to  some  thirty  diameters  in  the  experiment 
of  Chang)— the  short  inlet  tangent  resulted  to  a  developing  flow  at  the  beginning  of  the  bend  with 
the  thickness  of  the  boundary  layer  being  approximately  15  percent  of  the  hydraulic  diameter. 
Their  data  indicates  that,  despite  the  thinner  boundary  layer,  the  S-like  structure  of  the  streamwise 
velocity  orofiles  were  as  notable  in  their  measurements  as  in  those  of  Chang.  They  also  reported 
the  generation  of  a  very  strong  secondary  motion,  which,  by  135°  around  the  bend,  appeared  to 
have  broken  down  into  a  chaotic  pattern. 

Data  from  most  of  the  above  described  experiments  have  been  used,  over  the  last  ten  years, 
as  test  cases  for  validating  numerical  methods  and  turbulence  models.  Humphrey  et  al.  (1981) 
presented  calculations  for  their  experimental  configuration  using  the  standard  k-e  model  in 
conjunction  with  the  wall-functions  approach.  The  same  configuration  was  one  of  the  test  cases 
selected  for  the  1980-81  Stanford  Conference  on  Complex  Turbulent  Flows  (Kline  et  al.,  1981) 
and  was  computed  by  several  participants.  In  most  cases  good  agreement  between  computed  and 
measured  mean-velocity  components  was  obtained  for  the  first  half  of  the  bend  but  significant 
deviations  were  present  farther  downstream.  The  configuration  of  Taylor  et  al.  (1982)  has  been 
computed  by  Buggeln  et  al.  (1980)  (one-equation  turbulence  model),  Govindan  et  al.  (1991) 
(algebraic  turbulence  model)  and  Kunz  and  Lakshminarayana  (1991)  (k-e  model  with  wall- 
functions).  All  the  numerical  predictions  for  this  case  are  in  better  overall  agreement  with  the 
measurements  than  the  corresponding  results  for  the  configuration  of  Humphrey  et  al.,  despite  the 
fact  that  the  same  bend  was  used  for  both  cases.  This  is  mainly  due  to  the  thinner  boundary  layer 
present  at  the  inlet  of  the  bend  in  the  study  of  Taylor  et  al.  which  resulted  in  weaker  secondary 
motion  than  measured  by  Humphrey  et  al.  The  180°  bend  configuration  of  Chang  has  been 
studied  computationally  by  Chang  (1983),  Johnson  (1984),  Birch  (1984)  and,  more  recently,  by 


8 


Choi  et  al.  (1990).  Chang,  Johnson,  and  Birch  used  in  their  calculations  the  standard  two- 
equation  k-e  turbulence  model  in  conjunction  with  wall  functions;  Chang  also  carried  out 
calculations  using  an  algebraic  Reynolds  stress  model  but  he  was  unable  to  obtain  a  converged 
solution.  All  three  computations  failed  to  reproduce  significant  features  of  the  flow  field,  such  as 
the  S-shaped  profiles  of  the  mean  axial  velocity  component  caused  by  the  strong  secondary 
motion.  In  the  more  recent  study,  Choi  et  al.  showed  that  the  inability  of  the  computations  to 
reproduce  correctly  the  strength  of  the  secondary  flow  and  its  interaction  with  the  axial  flow  can  be 
attributed  primarily  to  the  inadequacy  of  the  wall-function  approach  for  complex,  three- 
dimensional,  non-equilibrium  shear  flows.  They  carried  out  calculations  all  the  way  to  the  wall  by 
replacing  the  wall  functions  with  a  fine  mesh  across  the  sublayer  and  te.sted  both  the  k-e  model  and 
an  algebraic  Reynolds- stress  (ASM)  model.  For  both  models,  the  near-wall  region  was  modeled 
by  the  van  Driest  mixing-length  formula.  Their  k-e  calculations  showed  significant  improvements, 
as  compared  to  previous  computations  (Chang,  Johnson,  Birch),  in  the  agreement  between 
measurements  and  computations.  Better  agreement  for  the  axial  velocity  profiles  and  the 
distribution  of  the  Reynolds  stresses  was  obtained,  however,  when  the  algebraic  Reynolds-stress 
model  was  employed.  Despite  marked  improvements  in  their  results,  significant  discrepancies 
between  experiment  and  calculations  still  remained.  However,  the  success  of  their  modelling 
refinement  sequence  indicates  that,  for  future  improvements  in  the  modelling  of  such  flows, 
research  efforts  should  focus  in  the  development  of  turbulence  models  which  account  for  non¬ 
isotropic  effects  and  at  the  same  time  provide  accurate  resolution  of  the  near-wall  layer. 

The  most  recent  experiment  listed  in  Table  2  is  that  of  Kim  (1991),  who  made  measurement 
of  pressure,  mean  velocity,  and  Reynolds  stresses  in  developing  boundary-layer  flow  in  a 
rectangular  duct,  of  aspect  ratio  six,  turning  through  an  angle  of  90'’.  There  are.  to  be  sure,  several 
previous  measurements  in  rectangular  ducts  of  different  aspect  ratio  but,  to  the  authors' 
knowledge,  most  of  these  were  concerned  with  .secondary  motion  (of  the  .second  kind)  in  fully- 
developed  flow  in  straight  ducts.  In  particular,  none  has  considered  developing  flow  in  a  curved 
duct.  Kim's  measurements  are  particularly  noteworthy  because  they  document  the  development  of 
the  pressure-driven  secondary  motion  in  boundary  layers,  and  the  formation  of  longitudinal 
vortices  within  the  boundary  layers  as  a  result  of  this  .secondary  motion.  His  data  also  reveal  the 
action  of  the  concave  and  convex  surface  curvatures  on  the  boundary  layer  velocity  and  turbulence 
profiles.  Kim  made  extensive  comparisons  between  measurements  and  calculations.  The 
calculations  were  performed  with  one  of  the  numerical  methods  used  in  the  present  work.  As  the 
method  will  be  de.scribed  in  this  report,  discussion  of  Kim's  results  is  postponed  until  a  later 
section. 

Measurements  in  turbulent  flow  in  curved  circular  ducts  appear  to  be  relatively  sparse  in 
spite  of  the  varied  practical  applications.  The  first  detailed  measurements  were  reported  by  Rowe 
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(1970)  who  used  Pitot  probes  to  measure  the  total  pressure,  mean  velocity,  and  yaw  angles  in 
several  cross-sections  of  a  180®  pipe  bend  of  mild  curvature.  His  data  indicate  that  the  crossflow 
reaches  a  maximum  at  a  bend  angle  of  30®  and  then  decreases  to  a  near  constant  value.  Rowe 
explained  this  phenomenon  by  examining  the  gradients  of  the  total  pressure  and  using  inviscid 
theories  to  predict  the  evolution  of  the  secondary  motion.  In  a  later  work,  Enayet  et  al.  (1982b) 
used  LDV  to  measure  the  streamwise  components  of  the  mean  and  fluctuating  velocity  in  a  90® 
strongly  curved  pipe  bend.  A  very  detailed  set  of  data  for  turbulent  pipe-bend  flow  was  reported 
by  Azzola  and  Humphrey  (1984)  who  measured  both  the  .streamwise  and  circumferential 
components  of  the  mean  and  fluctuating  velocity  in  a  180®  bend.  The  major  features  of  this  flow 
are  the  reversals  that  the  secondary  motion  undergoes  for  bend  angles  larger  than  90®,  a 
phenomenon  which  was  found  to  be  essentially  independent  of  the  Reynolds  number,  and  the  large 
levels  of  turbulence  anisotropy  arising  everywhere  in  the  bend  and  in  the  downstream  tangent. 
Comparison  of  this  pipe-bend  flow  with  the  corresponding  square-duct  flow  (Chang,  1983) 
reveals  that  the  magnitudes  of  the  .secondary  motion  are  in  general  lower  for  the  pipe  bend. 
Moreover,  there  are  no  crossflow  reversals  in  the  square  bend  case  which  demonstrates  that, 
despite  superficial  similarities,  the  two  flows  may  be  quite  different.  The  most  recent  set  of 
measurements  was  reported  by  Anwer  et  al.  (1989)  who  carried  out  mean  flow  and  turbulence 
measurements  through  a  180®  pipe  bend  and  in  the  downstream  tangent.  Their  measurements 
include  the  three  components  of  the  mean  velocity  and  the  six  Reynolds  stress  components. 

The  first  turbulent  flow  calculation  of  a  curved  pipe  flow  was  reported  by  Patankar  et  al. 
(1975),  who  presented  solutions  for  the  mildly  curved  bend  measured  by  Rowe.  Their  solutions, 
obtained  with  the  standard  k-e  model  with  wall  functions,  were  quite  successful  in  reproducing  the 
general  features  of  the  flow  observed  in  the  measurements.  However,  the  shapes  of  the  computed 
contours  of  constant  velocity  head,  which  are  not  as  distorted  as  the  corresponding  measured  ones, 
indicate  that  the  calculations  did  not  predict  accurately  the  development  and  decay  of  the  secondary 
motion.  It  should  be  pointed  out,  however,  that  the  data  of  Rowe  do  not  include  quantitative 
measurements  of  the  secondary  motion  and  thus  no  direct  comparisons  could  be  made.  In  a  more 
recent  study,  lacovides  and  Launder  (1984)  calculated  the  90®  pipe  bend  of  Enayet  et  al.  (1982b). 
They  employed  a  two-layer,  k-e  model  with  van  Driest’s  eddy-viscosity  formula  for  the  near-wall 
layer.  The  agreement  of  their  calculations  with  the  experimental  data  was  broadly  satisfactory  for 
the  streamwise  velocity  profiles,  but  the  boundary  layer  on  the  inside  of  the  bend  did  not  proceed 
as  far  toward  separation  as  indicated  by  the  measurements.  Moreover,  the  calculated  rate  of 
recovery  of  the  flow  in  the  downstream  tangent  was  slower  than  measured.  The  1 80®  pipe  bend  of 
Azzola  and  Humphrey  (1984)  was  computed  by  Azzola  et  al.  (1986),  using  the  same  method  as  in 
lacovides  and  Launder  (1984).  Their  calculations  for  the  development  of  the  streamwise  mean 
velocity  component  through  the  bend,  along  the  radius  at  90®  from  the  symmetry  plane,  were  in 
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good  agreement  with  the  measurements.  However,  significant  discrepancies  between  experiment 
and  calculations  were  present  in  the  corresponding  profiles  of  the  circumferential  mean  velocity 
component  between  bend  angles  of  45®  and  180®,  as  well  as  in  the  rate  of  recovery  of  the  flow  in 
the  downstream  tangent.  In  a  very  recent  study,  Lai  (1990)  computed  the  180®  pipe  bend  of 
Anwer  et  al.  (1989)  using  a  full  Reynolds-stress  transport  model  including  the  direct  modelling  of 
the  near  wall  flow.  His  calculations  emphasized  the  secondary  flow  patterns  in  curved-pipe  flows. 
Based  on  analytical  arguments  with  the  vorticity  transport  equation,  as  well  as  on  his  numerical 
calculations,  he  showed  that  a  turbulence-driven  secondary  motion  is  present  near  the  outer  bend  of 
the  curved  pipe. 

II. 4  Experiments  in  turbulent  flow  in  transition  ducts 

The  term  transition  duct  denotes  an  internal  flow  configuration  whose  cross-sectional  shape 
changes  in  the  streamwise  direction,  for  example,  from  circular  to  rectangular  (CR),  or  vice  versa. 
Transition  ducts  can  be  either  straight  or  curved  and  are  encountered  in  aircraft  propulsive  systems 
(as  inlet  and  exhaust  nozzles  of  jet  engines),  wind  tunnels  (wind  tunnel  contractions  and  diffusers), 
and  hydraulic  turbine  systems  (draft  tubes).  Turbulent  flow  measurements  through  transition  ducts 
have  been  reported  by  Mayer  (1939),  Taylor  et  al.  (1981),  Burley  and  Carlson  (1985),  Patrick  and 
Me  Cormick  (1987,  1988)  and,  more  recently,  by  Miau  et  al.  (1990),  and  Davis  and  Gessner 
(1992).  All  these  experiments,  which  have  focused  exclusively  on  straight  transition  ducts,  are 
summarized  in  Table  3. 

The  first  turbulent  flow  measurements  in  a  transition  duct  were  reported  by  Mayer  (1939) 
who  measured  the  flow  through  two  rectangular-to-circular  (and  vice  versa)  transition  ducts  of 
constant  cross-sectional  area.  His  measurements  included  streamwise  static  pressure  variations 
along  the  duct  walls,  total  pressure  contours  and  mean  velocity  components.  The  next  set  of 
measurements  was  reported  four  decades  later  by  Taylor  et  al.  (1981)  who  carried  out  mean 
velocity  measurements  through  a  square-to-circular  transition  duct  whose  cross-sectional  area 
reduced  by  21.5%  over  the  transition  length.  These  two  experimental  studies  demonstrated  the 
impact  of  the  transition  length  on  the  streamwise  flow  development  and  showed  that  a  moderate 
secondary  motion  (not  exceeding  10%  of  the  bulk  velocity)  can  induce  significant  distortion  of  the 
streamwise  flow.  In  a  more  recent  study,  Burley  and  Carlson  (1985)  carried  out  pressure  loss 
tests  for  transonic  flow  through  five  different  CR  transition  ducts  in  order  to  study  the  effect  of 
transition  length,  cross-sectional  shape  and  inlet  swirl  on  the  duct's  performance.  Their 
measurements  indicated  that  short  transition  lengths  (less  than  0.75  hydraulic  diameters)  can  induce 
large  regions  of  separated  flow.  They  also  found  that  the  inlet  swirl  can  enhance  the  performance 
of  low  pressure  ratio  transition  ducts. 
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The  first  turbulence  measurements  for  transition  duct  configurations  were  reported  by 
Patrick  and  McCormick  (1987,  1988).  They  carried  out  measurements  for  two  CR  transition 
ducts:  i)  the  first  with  an  aspect  ratio  of  three  and  transition  length  equal  to  one  diameter,  and  ii)  the 
second  with  an  aspect  ratio  of  six  and  transition  length  equal  to  three  diameters.  The  cross- 
sectional  area  ratio  remained  constant  for  the  first  duct  but  it  increased  to  1 . 1  at  the  midpoint  and 
then  decreased  back  to  1.0  for  the  second  duct.  Their  measurements  include  mean  velocity 
components  and  normal  Reynolds  stresses  at  the  inlet  and  exit  planes.  More  recently  Miau  et  al. 
(1990)  carried  out  turbulent  flow  measurements  through  three  CR  transition  ducts  of  constant 
cross-sectional  area.  They  reported  measurements  of  mean  velocity  components,  static  pressure 
distributions  and  normal  Reynolds  stress  distributions.  Miau  et  al.  used  their  data  to  evaluate  the 
different  terms  in  the  transport  equation  for  the  mean  streamwise  vorticity  component  and  showed 
that  the  generation  of  the  secondary  motion  is  due  to  the  transverse  pressure  gradients  induced  by 
the  rapid  geometrical  changes.  Perhaps  the  most  complete  set  of  measurements  for  a  transition 
duct  configuration  is  the  very  recent  work  of  Davis  and  Gessner  (1992)  who  measured  the  flow 
through  a  CR  transition  duct  (overall  length-to-diameter  ratio  of  4.5,  aspect  ratio  of  3.0  at  the  exit 
plane,  cross-sectional  area  ratio  increasing  to  1.15  at  midpoint  and  then  decreasing  back  to  1.0  at 
exit).  They  made  very  detailed  mean  flow  and  turbulence  measurements  (including  all  the  six 
components  of  the  Reynolds  stress  tensor)  at  several  cross-sections  within  the  transition  region  as 
well  as  at  the  exit  plane  of  the  duct.  Their  results  show  that  the  curvature  of  the  duct  walls  induces 
a  relatively  strong  pressure-driven  secondary  motion  which  significantly  distorts  both  the 
streamwise  mean-velocity  and  Reynolds-stress  fields.  In  addition,  careful  analysis  of  their  near¬ 
wall  data  indicated  that  the  extent  of  the  law-of-the-wall  behavior  is  diminished  in  regions  where 
streamwise  wall  curvature  effects  influence  the  flow  development. 

Finally,  we  note  that  data  on  turbulent  flow  in  curved  transition  ducts,  with  significant  area 
changes,  such  as  a  hydroturbine  draft  tubes,  are  quite  sparse,  and  certainly  neither  well 
documented  nor  detailed  enough  to  serve  the  purpose  of  evaluating  numerical  methods. 
Nevertheless,  we  shall  refer  to  the  available  data  in  a  later  section. 

III.  SCOPE  OF  THE  PRESENT  WORK 

The  foregoing  review  of  previous  research  on  flow  in  curved  ducts  clearly  indicates  that  the 
related  literature  is  vast  and  still  growing.  The  subject  continues  to  be  of  great  current  interest 
because  there  are  a  number  of  issues  that  remain  to  be  settled.  In  the  case  of  laminar  flow,  the 
central  concern  is  that  of  developing  numerical  methods  for  the  solution  of  the  Navier-Stokes 
equations  that  can  accurately  capture  the  various  phenomena  that  are  observed  in  the  experiments. 
There  exist  significant  differences  among  calculations  made  by  different  numerical  methods  with 
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identical  geometries  as  well  as  initial  and  boundary  conditions.  Such  differences  are  also  ob.served 
in  turbulent  flow.  In  this  case,  however,  experiments  reveal  a  much  more  complex  response  of  the 
flow  to  curvature,  and  the  numerical  problem  is  compounded  by  the  differences  in  turbulence 
models  that  are  used  for  closure  of  the  Reynolds-averaged  Navier-Stokes  equations.  The  initial 
and  boundary  conditions  employed  for  the  additional  turbulence-model  equations,  and  the  manner 
in  which  these  equations  are  solved,  influence  the  results.  In  both  laminar  and  turbulent  flows, 
there  are  a  number  of  other  problems,  such  as  grid  topology  and  grid  generation,  that  must  be 
resolved  when  existing  methods  are  extended  and  applied  to  complex  duct  geometries  of  practical 
interest.  Thus,  current  research  is  driven  by  the  need  for  accurate  numerical  methods,  on  the  one 
hand,  and  improved  turbulence  models,  on  the  other.  These  two  aspects  become  even  more  critical 
when  numerical  methods  are  u.sed  for  the  prediction  of  the  flow  in  curved  ducts  of  arbitrary  and 
changing  cross  section. 

The  present  work  repre.sents  the  first  pha.se  of  a  broad  re.search  effort  that  aims  to  develop 
numerical  methods  suitable  tor  predicting  the  flow  through  the  passages  of  hydraulic 
turbomachinery.  Passages  such  as  turbine  scrolls  and  draft  tubes  are  strongly  curved  ducts  of 
arbitrary  and  changing  cross-section.  In  view  of  the  inadequacy  of  state-of-the-art  numerical 
methods  and  turbulence  models  to  predict  flows  through  simply-curved  ducts  of  square  and 
circular  cross  section,  it  is  not  surprising  that  the  flow  in  ducts  with  the  geometric  complexity  of  a 
typical  hydroturbine  lies  beyond  present  capabilities.  Our  initial  objective  herein  is  to  assess  the 
performance  of  the  numerical  methods  and  turbulence  models  currently  available  at  the  Institute  in 
the  context  of  curved  duct  flows  and  identify  the  focal  areas  of  our  subsequent  research  efforts. 

We  first  present  the  Reynolds  averaged  Navier-Stokes  equations,  along  with  the  turbulence 
closure  equations,  in  Canesian  coordinates.  Subsequently,  we  de.scribe  and  examine  in  parallel 
two  numerical  methods,  namely,  the  one  developed  by  Patel  and  his  co-workers  (see  Kim.  1991), 
which  will  be  referred  to  as  Method  I,  and  the  other  developed  by  Sotiropoulos  ( 1991 ),  referred  to 
as  Method  II.  We  then  proceed  to  carry  out  a  series  of  laminar  and  turbulent  flow  calculations  in 
ducts  of  regular  (square,  rectangular,  and  circular)  cross-section  and  to  compare  the  results  with 
experimental  data.  The  laminar  calculations  offer  a  very  good  indication  of  the  spatial  resolution  of 
the  numerical  methods  while  the  turbulent  calculations  allow  the  evaluation  of  the  turbulence  model 
of  Chen  and  Patel  (1988)  which  is  incorporated  in  both  methods.  Finally,  we  report  results  of 
calculations  made  for  the  straight  transition  duct  of  Davis  and  Gessner  (1992),  and  a  draft  tube 
configuration  derived  from  the  one  at  the  Norris  Dam  in  Tennessee.  These  are  perhaps  the  first 
such  results  to  be  obtained  with  a  fine  grid,  resolving  the  flow  all  the  way  to  the  wall,  in  the  case 
of  the  draft  tube,  qualitative  comparisons  are  made  between  the  present  results  and  those  of  some 
recent  calculations  for  similar  draft  tubes  made  by  Vu  and  Shyy  (1990)  and  Agouzoul  et  al.  (1990), 
who  have  employed  rather  coarse  grids.  The  present  solutions  reveal  a  wealth  of  detail  that  is 
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difficult  to  capture  with  coarse-grid  calculations,  and  offer  the  opportunity  to  speculate  on  advances 
that  must  be  made  to  develop  a  predictive  capability  for  turbulent  flow  in  hydroturbine  ducts  and 
similar  geometries  encountered  in  numerous  other  applications. 

IV.  GOVERNING  EQUATIONS  IN  CARTESIAN  COORDINATES 

The  Reynolds-averaged  Navier-Stok.es  equations  for  an  incompressible  fluid,  non- 
dimensionalized  by  the  fluid  density  p,  reference  velocity  Uq,  and  reference  length  Lq,  in  Cartesian 
coordinates  xj  (=  xi,X2,X3)  are 


dt  ^  3xj  dxj  Re  3xj3xj  9xj  ^ 


and  the  conservation  of  mass  (continuity)  condition  is 


9xi 


=  0 


(1) 


(2) 


Here,  V,  is  the  mean  velocity  component  in  x\  direction,  vj  denotes  the  corresponding  fluctuating 
component,  an  overbar  denotes  an  ensemble  average,  and  Re  is  the  Reynolds  number  (Re  = 
UqLo/v  ),  where  v  is  molecular  viscosity.  The  Reynolds  stresses  are  related  to  the  mean  flow 
through  the  eddy  viscosity  Vj  by 


-  VjVj  =  Vt 


9xj  9x, 


(3) 


where  the  Kronecker  delta  function  is  given  by 


i  =j 

0,  i  5^=  j 


The  eddy  viscosity  Vt  is  related  to  the  turbulent  kinetic  energy  k  and  its  rate  of  dissipation  e  by  the 
Prandtl-Kolmogoroff  formula 


and  k  and  e  are  determined  by  the  modelled  transport  equations 


14 


|.+v,|!^=±[(x+i)|il  +G-£ 

ot  ^  dXj  dXj  [iRe  Gk/dxj 


^  +  V 
at  ^dxj 


^  a  r/ 1  ^  vt  \  ae 

axj  \Re  ag/axj 


.Q,£a.Q.L 


where  the  generation  of  turbulent  kinetic  energy,  G,  is 


^  _ av, 

G  =  -ViV,^=v, 


av,  av. 


[dxj  ax. 


dxi 


(5) 

(6) 

(7) 


The  constants  used  in  this  "standard"  k-e  model  are  =  0.09,  Cei  =  1.44,  Ce2  =  1-92,  =  1.0 

and  Gg  =1.3. 

In  the  two-layer  turbulence  model  of  Chen  and  Patel  (1988),  which  is  employed  in  both 
methods  here,  the  flow  domain  is  divided  into  two  regions  -  the  inner  layer  and  the  outer  layer. 
The  inner  layer  includes  the  sublayer,  the  buffer  layer  and  a  part  of  the  fully-turbulent  or 
logarithmic  layer.  The  simple  one-equation  model  of  Wolfshtein  (1969)  is  modified  and  employed 
to  account  for  the  wall-proximity  effects,  whereas  the  standard,  two-equation  k-e  model  is  used  in 
outer  layer. 

The  one-equation  near-wall  model  requires  the  solution  of  only  the  turbulent  kinetic  energy 
equation  (5)  in  inner  layer.  The  rate  of  energy  dissipation  in  this  region  is  specified  by 


and  the  eddy  viscosity  is  obtained  from 

v,=  C^VF('^  (9) 

where  the  length  scales  4  ^nd  fp  contain  the  damping  effects  in  the  near-wall  region  in  temis  of 

Vic  y 

the  turbulence  Reynolds  number  Ry  (= - : 

^  V 


V  y 


1  -  exp 


(10) 


4  =  C,y 


1  -  exp 


(11) 
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Note  that  both  and  4  become  linear,  and  approach  y  with  increasing  distance  from  the  wall. 
The  constant  is  given  by 

C,  =  kC^-3/4  (12) 


K  (=0.418)  being  the  von  Karman  constant,  to  ensure  a  smooth  eddy-viscosity  distribution  at  the 
junction  of  the  inner  and  outer  layer.  In  addition,  Ag  =  2C^  is  assigned  so  as  to  recover  the  proper 

asymptotic  behavior 


e  = 


(13) 


in  the  sublayer.  The  third  parameter  =  70  was  obtained  from  numerical  tests  to  recover  the 
additive  constant  B  =  5.45  in  the  logarithmic  law  in  the  case  of  a  flat-plate  boundary  layer.  In  the 
outer  layer,  beyond  the  near-wall  layer,  the  standard  k-e  model  is  employed  to  calculate  the  velocity 
field  as  well  as  the  eddy  viscosity. 

V.  DESCRIPTION  OF  METHOD  I 


The  numerical  method  developed  by  Patel  and  his  co-workers  (see,  for  instance  Kim, 
1991),  which  for  convenience  will  be  referred  to  as  Method  I,  is  described  first.  Method  1  utilizes 
the/n//  transformation  approach,  which  transforms  both  dependent  (velocity  and  Reynolds  stress 
components)  and  independent  (spatial  coordinates)  variables  to  generalized  curvilinear 
components.  In  the  following  sections,  the  governing  equations  in  generalized  curvilinear 
coordinates,  expressed  in  terms  of  the  physical  contravariant  velocity  components,  are  presented 
and  the  overall  solution  procedure  is  described. 

V.l  Governing  equations  in  general  curvilinear  coordinates 

The  transformation  of  the  equations  from  Cartesian  coordinates,  to  the  general 
nonorthogonal  coordinates,  is  facilitated  by  the  use  of  the  metric  tensor  in  the  coordinates 


dx*- 


(14) 
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and  its  inverse  matrix  g‘J  given  by 


8^^  ~  g  (SktnSln  ■  8kn8lm) 


(15) 


where  the  indices,  grouped  as  (i,k,l)  and  (j,m,n),  follow  a  cyclic  order.  The  Jacobian  of  the 
transformation  is  the  square  root  of  the  determinant  of  the  gjj  matrix,  j2  =  Igjjl  =  g,  and  the 
Christoffel  symbols  of  the  second  kind  are  related  to  the  metric  coefficients  by 


rlnn  ~  2 


/^8jm  ^  ^§nj 

\  a^ 


(16) 


To  carry  out  the  transformations,  the  partial  derivatives  are  first  changed  into  covariant 
derivatives.  Next,  the  vector  components  of  velocity  are  replaced  by  the  contravariant 
components,  and  the  indices  of  covariant  terms,  like  the  gradient  of  a  scalar,  are  rai.sed  by 
multiplying  by  the  geometric  coefficient  g*).  The  resulting  equations,  in  general  curvilinear 
coordinates,  are: 


continuity  equation 
V%=0 

momentum  equation 

^  +  (V”  -  g”"  V,  „|V' -  v,,„g™ V",,„ 

=  .g,n,|£^ai^qX.v,)g»n(v.4 


where  use  is  made  of  the  eddy-viscosity  relation 


Igir^ni  ^  gms^i  _  2k_  gim 


k-e  model  equations 


(17) 


(18) 


(19) 


(20) 

(21) 
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in  which 


G  =  -  gij  v‘v-  VJ,^  =  Vt(v-,V*.^  +  gyg-"V,„,VJ,n) 


(22) 


In  the  above,  the  covariant  derivative  of  a  scalar  is  defined  by 
9k 


km  —  ' 


9^" 


(23) 


the  covariant  derivative  of  a  contravariant  component  of  a  vector  is  given  by 

v‘  =-^  +  rLuV‘^ 


9^" 


(24) 


and  the  covariant  derivative  of  a  covariant  component  of  a  vector  is  given  by 
9(k,in)  r-s  /I  \  9  /  9k  \  9k 


/I,  \  o  \ 

...n  *  1  mn(K.s)  -  I  -  1  n 


9^" 


mn 


9^"\9^'"/  9^‘ 


(25) 


The  covariant  derivative  of  a  mixed  second-order  tensor  is  given  by 


(V\„,),n  =  +  ^"ns(V^^,)  -  r1nn(vM 


(26) 


which  results  in 


g'""(V',n,).n  =  vV  +  2C-^+V* 


_L  rSmr"  _i_  /~'imr^ 
~  ^  '-'S  •  mn  ^  '^n  ‘  m: 

[9^'" 


(27) 


where  =  g"’"rj,5  >  and  is  the  Laplacian  of  a  scalar,  defined  by 


=  gnin _ ^ _ +  fni_^ 

9^m^n 


(28) 


and 


f^  =  -g-r?^  =  |-^(Jg^'") 
J  9^^ 


(29) 
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The  contravariant  components  of  velocity  vectors  do  not,  in  general,  have  measurable 
dimensions.  Therefore,  following  Trusdell  (1953),  the  so-called  physical  components  of  the 
contravariant  components,  defined  by 


V(i)  =  V^V 


(30) 


are  introduced.  Then,  the  final  equations,  which  are  solved,  are  obtained  by  substituting  (30)  into 
equations  (17)  through  (21).  Thus,  we  have  the 


continuity  equation 


=  0 


(31) 


momentum  equarion 

av(i) 


at 


■  +  gm/n 


V(m)^  +  D^(i)V(m)V(i)  +  g|/2g>{?rLV(m)V(n) 


-  ^  [gmn^  +  g-n‘,/2g|/2g.m^  +  E"’(i)V(i)  +  g-„1^2y;.2E.(n)V(n 


a^"  a^' 


a^ 


+  gj/Vs^Q"  +  Q‘)V(S)1  + 


-I- 


?V(i)  +  2  E"'(i)^  +  P(i)V;il  +  sJ/2s:,'^2l2  +  Qi,V(n)\ 

ar  I  ar  I. 


(no  sum  on  i) 


(32) 


k-c  moH^  ^uadons 


ak 

at 


+  girl^^ 


V(m)  -  —  g'"” 


avt 


cJk  3^" 


ak 


=  (-L-h  A|v\-HG-e 

\Re  aj 


V(m)  -  ^  g’ 


jmn 


av, 


arj3^ 


ae 


iRe  CTe/  k  k 


(33) 


(34) 


where 

G  =  -  gij  V‘V"1  VJ.m  =  V,  (V'",,V',m  +  gijg"’"  V'.niVJ,n) 
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with 


av(i) 


V‘ n,  =  zP  ^  D„,(i)V(i)  +  g-,i/2  rLv(s) 


3^ 


and  other  geometry-related  terms  are 


Dn,(i)  =  - 


1  3gii 

2gi. 


=  gmn 

C^"  =  g"’"  Cs 

P(i)  =  + [D^(i)  -H  r"J  E-(i) 

dC 

ar*"’  r  cl  c 

Qn  =  ^  +  [2DUi)  +  rij  cr  +  Tinker 

dc 


(35) 


V.2  Discretization  of  the  transport  equations 

The  three-dimensional  form  of  the  transport  equations  of  momentum  and  turbulence 
properties  are  discretized  using  the  finite-analytic  method  developed  by  Chen  and  Chen  (1984)  and 
modified  by  Patel,  Chen  and  Ju  (1988).  The  dimensionless  forms  of  the  momentum  and  k-e 
equations  are  arranged  into  the  form  of  a  general  convection-diffusion  equation  as  follows: 

+  g22(t)^^  +g33(t>^i;  =  2C<„(t>4  +  2B<,(t)T,  2A^(t)i;  -t-R*(|),  + 

(4d) 

where  the  subscript  indicates  a  partial  derivative  with  respective  to  time  or  one  of  the  general 
coordinates  (^*,^2  ^3)  =  (^,r|,Q.  The  variable  <{)  represents  one  of  the  following  quantities:  [V(l), 
V(2),  V(3),  k,  £]. 

When  (j)  represents  one  of  the  velocity  components,  [V(1),V(2),V(3)]  =  (U,V,W),  say,  the 
coefficients  and  source  terms  derived  from  the  momentum  equations  are 


2C0  -  Rrf 


8i'pV(l)-gi"^ 

3^ 


g 


ii 


9vt 


-[f*  -h2Ei(i)  +  2Ci‘] 


(37) 
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2B^  =  r g-.HMl)  -  g-"  ^  -  g‘2  ^1  -  [f2  +  2E^(i)  +  2Cp] 

2A^  =  Rjgs'f  V(3)  -  g3"  ^  -  g‘3  ^1  -  [f3  +  2EHi)  +  2Cp] 

=  Rv'  =  Rv-  =  Rv^  =  -r-^ - 

Sj  =  -  NO(V(i)]  +  R,  gl/jfg™  -^p  + 

. 

T"  3^'  s^'"  °  ae  J 

-  R<,  -^[E'"(i)V(i)  +  gi^2gWE'(mlV(m)  +  g-nW2gl/2(c»  +  C!r'|V(n)  ] 

-  P(i)V(i)  +  g!«gs!iSQ;.,v(m)  -f  gingjiQac™^  +  gWgOi^njc™?^^ 
+  R,i,ghli<?[Dni(i)V(i)V(m)  +  g|/’ghy“rJnnV{m)V(n)] 


where  (i,j,k)  are  in  cyclic  order  and  the  diffusion  terms  arising  due  to  the  nonorthogonal 
coordinates  are  given  by  the  operator 

NO((!))  =  2[g>-(!)^^  +  g‘3(t)^^  +  g23(i,^  J 


The  various  geometric  coefficients  were  defined  in  the  previous  section.  Note  that  they  involve 
sums  over  the  indices  m  and  n.  The  coefficient  required  for  each  component  of  the  momentum 
equation  can  be  found  using 

(|,  =  V(1)  =  U  (i=l,j=2,k=3) 

(1)  =  V(2)  =  V  (i=2,j=3,k=l) 

,|)  =  V(3)  =  W  fi=3,j=l,k=2) 

The  coefficients  for  (()  =  k  or  e  are  derived  from  their  corresponding  non-dimensional 
equations  and  arranging  them  in  the  form  of  equations  (37)  to  (39)  and  (41).  This  yields 


2 


(43) 


2Ct  =  Rjg-,Y-V(l)-;Lgl.i^l.fl 

[  3|  J 


2B,  =  R  Jg2Y=V(2)  ^  g2"  ^1  -  fi 

I  3^  J 


2A,  =  Rjg3Y=V(3)-^g3"|^]-f3 

'-'0  dq 


Sk  =  -  NO(k)  -  Rk  (G  -  e) 


Sg  —  N0(£)  -  Rg  ^  (CgiG  -  Cg26) 


where  R<|)  is  the  effective  Reynolds  number  defined  by 


-L  =  -L  +  :^ 
Rk  Re  Ok 


Rg  Rc  Qig 


Now,  introducing  the  coordinate  stretching 
I  =  ^  T1*  =  ^  ^ 

Vip  Vy  (50) 

equation  (36)  reduces  to  the  standard  three-dimensional,  convective-diffusive  transport  equation 
described  in  Chen  and  Chen  (1984),  i.e., 

0^*^*  +  ~  2C(j)^*  +  +  2A({)(;*  +R(|)t  +  (S,j,)p 
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(52) 


C  _  (Co)P  ^  B  _  (B0)P  .  ^  _  (A0)P  ^ 

VgJ^ 


R  =  (R0)p 


for  a  numerical  element  with  dimensions 


ATI* 


AC 


(53) 


shown  in  figure  2,  where  subscript  P  denotes  a  quantity  evaluated  at  a  point  (node)  P. 

In  applying  the  finite-analytic  method  of  Chen  and  Chen  (1984)  to  the  linearized  equation 
(51),  the  most  general  version  for  a  three-dimensional  element  results  in  a  28-point  discretization 
formula.  While  such  a  scheme  will,  in  principle,  yield  better  results,  for  a  wide  class  of  problems 
it  suffices  to  use  a  simplified  method  to  save  computational  effort.  Here  the  hybrid  method,  which 
combines  a  two-dimensional  local  solution  in  the  TiC"plane  with  a  one-dimensional  solution  in  the 
^-direction,  is  employed.  A  brief  description  is  given  below. 

In  the  hybrid  scheme,  equation  (5 1 )  is  decomposed  into  one-  and  two-dimensional  partial 
differential  equations  as  follows 


2C(t)^*  -  -t-  R(t)t  +  (S^jp  =  G(C*,r|*,C*,t) 

-  2A(t)^*  =  G(C*,'n*,C*4) 


(54) 

(55) 


If  G  is  assumed  to  be  constant  in  each  local  element,  and  if  the  time  derivatives  are  approximated 
by  backward  finite-difference  formulas,  equations  (54)  and  (55)  reduce,  respectively,  to  the 
standard  one-  and  two-dimensional  convection-diffusion  equations  described  in  Chen  and  Chen 
(1984). 

The  general  solution  of  the  one-dimensional  equation  (54)  can  be  readily  obtained  as 
(|)  =  a(e2C^*  -  l)-i-b^*  +  c  (56) 

where  a,  b,  c  are  constants  to  be  determined  from  the  boundary  conditions.  By  substituting 
equation  (56)  into  (54),  the  source  function  G  at  the  center  node  P  is  obtained  as 


Gp  =  g  =  G(0,0,0,0)  =  (2C(j)^*  -  (})^*^*  +  R0I  +  S^)p 
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(57) 


=  (Cu  +  Cd)  ())?  -  Cu^u  -  Cd<1)d  +  —  (<|>p  -  0p  ^)  +  (S0)p 

Ax 


with 


=  - ,  C^  =  — -  (58) 

^  (’  sinh  Cf  °  Z'  sinh  C/ 

where  the  subscripts  P,  U,  and  D  denote  the  center,  upstream  and  downstream  nodes, 
respectively,  as  shown  in  figure  2.  The  superscript  (n- 1 )  refers  to  the  value  at  the  previous  time 
step,  and  Ax  is  the  time  step. 

By  specifying  boundary  conditions  for  each  element,  equation  (55)  can  be  solved 
analytically  by  the  method  of  separation  of  variables.  The  boundary  conditions  on  all  four 
boundaries,  r|*  =  ±  k  and  ^*  =  ±  h  of  the  transverse  section  (Ti(^-plane)  of  each  local  element  are 
specified  as  a  combination  of  exponential  and  linear  functions,  which  are  the  natural  .solutions  of 
the  governing  equation,  i.e., 

(t)(k,C*)  =  an  (e2AC*  -  l)  +  bnC*  +  Cn 
(t)(-k,;*)  =  a,(e2A;*-  i)  +  b,(;*  +  c, 

(()(ri*,h)  =  ae  (e^Bn*  -  l)  +  beH*  +  Ce 

(t)(tl*,-h)  =  aw  (e2BTi*  -  1 )  +  bwTl*  +  Cw 


where  a,  b,  c  are  again  constants  determined  by  the  nodal  values  along  the  boundaries  of  the 
element.  When  the  local  solution  thus  derived  is  evaluated  at  the  center  node  P  of  the  element,  the 
following  nine-point  discretization  formula  is  obtained: 

$p  =  Cne^ne  +  Cnw<(>nw  +  Cse4>se  +  Csw<l)sw 

+  CEC<t>EC  +  Cwc<f>wc  +  CNC<t>NC  +  Csc<}>SC  +  Cp  g  (59\ 


where 


Csc  = 


>Bk 


2  cosh  Bk 


Pa 


Cnc  =  e'2Bk  Csc 
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Cwc  = 


2  cosh  Ah 


Pb 


CeC  =  e'2Ah  Cwc 


Csw  = 


gBk  +  Ah 

4  cosh  Bk  cosh  Ah 


(I-Pa-Pb) 


CsE  =  Csw 

Cnw  =  Csw 

Cne  =  Csw 

Cp  =  1  .  P^)  =  kiupc  ,  1 .  Pj, 


with 


Pa  =  4E2  Ah  cosh  Ah  cosh  Bk  coth  Ah 


Pb  =  1  + 


Bh  coth  Bk  /p  |\ 
Ak  coth  Ah'  ^  ’ 


and 


E2=X 


-i-irKh 


"•■"1  [(Ah)2  +  ^  cosh  Va^  +  B^  +  k 


Finally,  by  substituting  the  inhomogeneous  term  g  from  equation  (57)  into  equation  (59).  the 
following  discretization  formula  is  obtained  for  an  unsteady,  three-dimensional  flow: 
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(  Cne<|>NE  +  Cnw(|>NW  CsEft>SE  +  Csw<|)SW 


1  +  Cp  (Cu  +  Cd  +  — ) 

'  Ax’ 

+  Cec^ec  +  Cwc<l)wc  +  Cnc<1>nc  +  Csr<|)sc 

+  Cp  |Cu<{)u  +  Cd<1)d  +  —  <t>P  * j  +  Cp  (S^)p } 

'  Ax  I  (60) 


or 


(j)p  = 


_ ! _ 

1  +  Cp  |Cu  +  Cd  +  — I 
'  Ax’ 


8 

{  Cnb0nb 
nb=l 


+  Cp  (  Cu(t)u  +  Cd(1)d  +  —  0P  '  I  +  Cp  (S^))p } 

Ax 


(61) 


where  the  subscript  nb  denotes  "neighboring"  nodes  (NE,  NW,  etc.).  It  is  seen  that  (})p  depends 
on  all  eight  neighboring  nodal  values  in  the  transverse  plane  as  well  as  the  values  at  the  upstream 
and  downstream  nodes  (<t)u,<|)D)  and  the  previous  time  step  .  When  the  cell  Reynolds 
number  2C  becomes  large,  Cy  approaches  ICK  =  2{C^)p  and  Cq  goes  to  zero  so  that  the 
coefficients  of  discretized  equations  turn  windward.  Similarly,  all  finite-analytic  coefficients  of 
transport  equations  automatically  control  the  direction  of  difference  and  weighting  of  neighboring 
nodes. 

As  equations  (61)  are  implicit  in  time  and  space,  their  assembly  for  all  elements  in  the 
solution  domain  results  in  a  system  of  simultaneous  algebraic  equations.  These  equations  are 
solved  by  the  tridiagonal-matrix  algorithm  in  conjunction  with  an  Altemating-Direction-lmplicit 
(ADI)  scheme.  For  steady-flow  calculations,  where  it  is  not  necessary  to  obtain  a  fully  converged 
solution  at  intermediate  stages,  only  ten  internal  iterations  are  used  during  each  time  step. 

V.3  Continuity  equation  and  pressure-velocity  coupling 

If  the  pressure  is  known,  equation  (61)  can  be  employed  to  solve  equations  (32)  for  three 
velocities,  and  (33)  and  (34)  for  k  and  e.  However,  the  pressure  is  not  known  a  priori  and  must  be 
determined  by  requiring  the  velocity  field  to  satisfy  the  continuity  equation  (31).  In  the  method  of 
Chen  and  Patel  (1989),  the  pressure  equation  is  derived  by  introducing  pseudo-velocities  at 
staggered  locations  while  maintaining  the  regular  grid  arrangement  for  all  the  transport  quantities. 
Figure  2  shows  the  locations  of  nodes  in  the  regular  grid  in  the  plane. 
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In  the  regular  grid  arrangement,  the  twelve-point  discretization  formula  (61)  for  the 
momentum  equations  gives  the  contravariant  physical  velocity  component  V(i)  at  P  as 


[V(i)]p  = 


1 


1  -I- Cp  (Cu  +  Cd  +  — 1  nb=l 
\  Ax/ 


{  S  CnbV(i)„b 


+  Cp  {CuV(i)u  +  CDV(i)D  +  — V(i)p  *}  +  Cp(Sv(i))p} 

Ax 


(62) 


Note  that  these  equations  contain  the  pressure  gradient  terms  inside  the  source  terms;  see  equations 
(41). 

The  actual  velocity  field  V(i)  in  the  momentum  equation  (32)  is  decomposed  into  a  modified 

pseudo-velocity  ^(i)  plus  the  pressure  gradient  terms  acting  in  the  direction  of  the  velocity 
component,  i.e.. 


V(i)  =  V(i)  - 


Cp  R  g!/^ 

1  -bCpjCu  +  CD+^j 


(63) 


V(i)  =  V(i)  -  E“  —  (no  sum  on  i) 

(64) 

The  continuity  equation  (3 1 )  yields 
[jgiY'v(i)]d-[jgiY'v(i)]u 
+  [jg2Y^V(2)]e-[jg-2Y2V(2)]w 

+  [jg3Y^V(3)]n-[jg3Y^V(3)]s  =  0  (65) 

Substitution  of  (64)  into  (65)  gives 

{[jgiY^E“]d  -b  [jgiY^E/  ‘L  -H  [jg2Y^E22],  +  [jg2Y2E22], 

+  [jg3V'E33]„+[jg3Y2E33]jpp 

=  [jgiY^E*  ‘]d  PD  +  [jg'iY^E  ‘  ‘]u  PU  +  [jg2Y^E22]g  PEC 

+  [jg2Y^E22]^  +  [jg^y2E33]^  pj^^+  [jg3y2E33]^  p^^  .  5 
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where 


D = [jgiY-v(i  )L + [jg-,V-v(i)L  +  [jgjV^vcaA 

+ [jgiY^wJw  +  [jgjY^voi  +  [jg-jY^vo)! 

With  a  regular  grid,  the  modified  pseudo-velocities  at  the  staggered  locations  cannot  be  obtained 
directly.  Therefore,  a  simple  one-dimensional  linear  interpolation  is  used  to  obtain  the  modified 
pseudo-velocities  at  the  required  staggered  nodes,  i.e., 

[v(l)]d  =  ^l[v(l)]D  +  [v(l)]pi 

[v(2)L  =  lj[v(2)]Ec+[v(2)]p) 

I  (OX) 

[v(3)L  =  1![v(3)]nc  +  [v(3)1,.) 

Z  PtC 


so  that  D  of  equation  (67)  can  be  expressed  in  terms  of  the  pseudovelocities  at  regular  nodes  as 
follows: 

D  =  ^  1[v(1)]d  -  [V(l)]u  +[v(2)]ec  -  [V(2)]wc  +  [v(3)]nc  -  [v(3)]sc) 

The  coefficients  (Jg'j  (Jg’V.,E22)g,  etc.  in  equation  (66)  are  also  calculated  by  one¬ 

dimensional  interpolation  in  the  same  way  as  in  equation  (68). 

The  solution  of  the  coupled  momentum  and  continuity  equations  involves  a  global  iteration 
process,  in  which  the  velocity-pressure  coupling  is  effected  by  predictor-corrector  steps.  In  the 
predictor  step,  the  pressure  field  at  the  previous  time  step  is  u.sed  in  the  solution  of  the  implicit 
momentum  equations  to  obtain  the  corresponding  velocity  field.  Since  this  velocity  field  does  not 
satisfy  mass  conservation,  a  corrector  step  is  needed.  In  the  corrector  step,  the  explicit  momentum 
equations  and  the  implicit  pressure  equation  are  solved  iteratively  to  ensure  that  the  continuity 
equation  is  satisfied.  This  procedure  is  almost  the  same  as  the  PISO  algorithm  of  Issa  (1985) 
except  for  some  minor  details.  In  the  PISO  algorithm,  the  pressure  equation  in  an  incremental 
form  is  solved  in  two  correction  steps.  The  present  procedure,  however,  solves  the  absolute 
(nonincremental)  pressure  equation  in  several  correction  steps,  the  number  of  which  is  specified  a 
priori.  In  the  present  algorithm,  the  finite-analytic  coefficients  are  not  updated  during  the  corrector 
steps,  while  in  the  PISO  algorithm  the  corresponding  convection-diffusion  coefficients  based  on 
finite-difference  formulation  are  updated  in  the  corrector  steps. 
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V. 4  Solution  procedure 

The  system  of  algebraic  equations  formed  by  (61)  for  transport  quantities  and  (66)  for 
pressure  is  solved  by  the  tridiagonal  matrix  algorithm  in  conjunction  with  the  ADI  (Alternating- 
Direction-Implicit)  scheme.  For  transient  problems,  where  the  initial  and  boundary  conditions  are 
properly  specified,  the  overall  numerical  solution  procedure  may  be  summarized  as  follows. 

1 .  Read  in  the  grid  and  calculate  the  geometric  coefficients  and  related  terms. 

2.  Specify  the  initial  conditions  for  velocity  and  turbulence  fields. 

3.  Give  an  initial  guess  for  pressure. 

4.  Calculate  the  finite-analytic  coefficients  for  the  transport  equations. 

5.  Solve  equations  (33)  and  (34)  for  turbulence  parameters  by  tridiagonal  matrix 
algorithm.  Update  the  eddy-viscosity  using  the  two-layer  turbulence  model. 

6.  Solve  the  momentum  equations  (32)  using  the  previous  pressure  field.  This  step 
constitutes  the  predictor  step  for  velocity  field. 

7.  Calculate  modified  pseudovelocities  from  equation  (64)  and  solve  the  pressure 
equation  (66). 

8.  Using  the  newly  obtained  pressure,  calculate  the  new  velocity  field  explicitly  from 
equation  (64).  This  completes  the  corrector  step. 

9.  Repeat  steps  7  and  8  a  spf’Citled  number  of  times  to  obtain  a  converged  solution. 

10.  Return  to  step  4  for  next  time  step. 

VI.  DESCRIPTION  OF  METHOD  II 

Method  il  is  based  on  the  one  developed  by  Sotiropoulos  (1991).  However,  during  the 
course  of  this  work,  a  number  of  new  features  and  capabilities  have  been  added,  including  an 
implicit  solution  procedure  for  the  pressure  equation  and  the  two-layer  turbulence  model  of  Chen 
and  Patel  (1988).  As  this  method  utilizes  the  partial  transformation  approach,  with  the  velocity 
components  expressed  in  Cartesian  coordinates,  the  appropriate  forms  of  the  Reynolds-averaged 
Navier-Stokes  equations,  and  the  equations  for  the  turbulence  model,  in  generalized  curvilinear 
coordinates  are  presented  once  again.  The  spatial  and  temporal  discretization  of  the  governing 
equations,  pressure-velocity  coupling  algorithm,  and  convergence  acceleration  techniques  are  then 
described  in  turn. 

VI. l  Governing  equations 

The  Reynolds-averaged  Navier-,Stokes  equations  (1)  and  the  equation  of  continuity  (2),  are 
transformed  to  generalized  curvilinear  coordinates  by  invoking  the  partial  transformation,  i.e.,  xi 
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— >  but  leaving  the  velocity  components  Vj  in  Cartesian  coordinates,  but  now  denoted  V'.  The 

transformed  governing  equations,  of  continuity  and  momentum,  read  as  follows: 


continuity  equation 


1  3(JV‘)  ^0 

momentum  equation 


(69) 


*  ‘  34' 


(70) 


In  the  above  equation,  Q  is  the  velocity  vector  Q  =  (  Vi,V2,V3  )^.  The  matrices  Aj  in  equation 
(70),  are  diagonal  matrices  defined  as  follows 


Aj  =  diag(VJ,VJ,VJ) 


(71) 


The  viscous  flux  vectors  E^j  which  appear  in  the  momentum  equation  (  70)  are 

E„,  =  (El„  E-;,  Ejf 


where 


Slj  —  ^X3^31 

S2j  =  ^i,Rl2+  ^^3^31 
S3j  =  ^i,Rl3+  ^X:R23 


k 


(72) 


The  source  vector  H  in  the  momentum  equation  (70)  is  given  by 
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„  3P  pk  3P  j.k  3P  -k 

34“  3^'’  3^“ 


The  equations  for  the  turbulent  quantities  k  and  £  are  written  in  generalized  curvilinear  as 


follows: 


3k  3k  1  3  f/  1  V,  \  3k  1  ^  „ 

— +  VJ—  F-  +  —  Jg‘^—  -G  +  e  =  0 


The  production  term  G  can  be  expressed  in  generalized  curvilinear  coordinates  as 


G=^v,  (R,j  +  Ri02 


where  repeated  indices  imply  summation  and  Ry  is  defined  in  equation  (72). 

VI.2  Spatial  discretization  of  continuity  and  momentum  equations 

The  momentum  equations  (70)  are  discretized  in  space,  on  a  non-staggered  mesh,  using 
three-point  central  finite  differencing  for  the  pressure  gradient  and  viscous  terms,  and  second  order 
upwind  finite  differencing  for  the  convective  terms.  The  upwind  differencing  of  the  convective 
terms  eliminates  the  need  for  adding  artificial  dissipation  terms,  to  the  right  hand  side  of  the 
momentum  equations,  to  stabilize  the  numerical  algorithm.  This  is  due  to  the  fact  that  a  fixed 
amount  of  dissipation  is  inherent  in  the  upwind  differencing. 

Referring  to  the  computational  cell  of  figure  3,  discrete  approximations  of  convective, 
pressure  gradient  and  viscous  terms  in  equations  (70)  are  as  follows: 


J.l,k 


3Xi  r  ~  (^xi)i,j,k  5^'  Pi.j.k 

a^L.k 
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where 


'_a_ 


V^  =  i[V±|v‘|],  (fori  =1,2,3) 

5?‘()i.j,k  =  ± — ^  [-3()i.j,k  +  4()i±i,jje  -  ()i±2.j,k] 

2A^ 

5^‘()i,j,k  =  ~[()i+l.j.k  '  ()i-l,j,k]  (XO) 

S^'Oi.j.k  =  — [()i+l/2,j.k  -  ()i-l/2.j.k] 

Similar  expressions  are  used  to  discretize  derivatives  with  respect  the  other  two  spatial  directions. 
In  all  the  above  equations,  the  metrics  and  the  Jacobian  of  the  geometric  transformation  are 
computed  at  the  (i,j,k)  nodes  using  three-point  central  differences.  To  compute  the  metrics  and  the 
Jacobian  at  the  half  nodes  ((i+l/2,j,k)  and  etc.),  where  they  are  needed  for  the  discretization  of  the 
viscous  terms,  a  simple  averaging  procedure  is  employed. 

The  continuity  equation  (69)  is  discretized  using  three-point  central  differences.  For 
convenience  we  define  the  discrete  divergence  operator  as  follows: 

=  (81) 

where  Q  is  the  cartesian  velocity  vector  and  V‘1  are  the  contravariant  components  of  the  velocity. 

VI.3  Temporal  discretization  of  continuity  and  momentum  equations 

The  system  of  the  discrete  continuity  and  momentum  equations  is  integrated  in  time  using 
the  four-stage,  explicit  Runge-Kutta  scheme.  Although  explicit  in  time,  the  Runge-Kutta  scheme- 
first  used  by  Jameson  (1981)  to  solve  the  compressible  Euler  equations— is  known  to  have  very 
good  error  damping  properties,  which  can  be  further  enhanced  by  employing  convergence 
acceleration  techniques  as  discussed  later  on.  Particularly  in  complex  three-dimensional  flow 
applications,  the  Runge-Kutta  scheme  can  be  very  competitive  to  approximate-factorization 
techniques,  since  it  is  fully  vectorizable  and  can  easily  take  advantage  of  parallel  processing 
capabilities  of  today's  supercomputers.  The  Runge-Kutta  scheme  has  been  adopted  by  several 
researchers  and  applied,  with  a  great  deal  of  success,  to  solve  compressible  viscous  and  inviscid 
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flows  (Jameson  and  Baker  (1983),  Chima  and  Yokota  (1990)).  Merkle  et  al.  (1986)  used  the  four- 
stage  Runge-Kutta  scheme  in  conjunction  with  the  artificial  compressibility  method  to  solve  the 
incompressible  Euler  equations,  while  Sotiropoulos  (1991)  applied  it  in  conjunction  with  the 
pressure-Poisson  approach  to  calculate  incompressible  turbulent  flows  past  ship  hulls. 

Following  Sotiropoulos  (1991),  the  Runge-Kutta  scheme  is  applied  to  the  system  of  the 
governing  equations  (69)  and  (70)  as  follows  (for  =  1, 2,3,4): 

DIV[Q^_^J=0  (82) 


Q  ‘  .  =  Q"  .  -  a,  At  .  RHS^:‘ 


(83) 


In  the  above  equations,  the  superscript  "n"  denotes  the  time  step  at  which  the  .solution  is  known, 

while  the  superscript "  denotes  an  intermediate  time  level  (or  iteration  level)  used  to  advance  the 

solution  from  time  step  "n"  to  time  .step  "n+l"  (we  designate  =  Q"  for  t=0  and  Q<  =  Q"’'’'  for 
^=4).  For  the  four-stage  .scheme,  the  coefficients  a/s  are:  1/4,  1/3,  1/2  and  1  for  <  =  1.2. 3.4.  in 

sequence.  The  RHS  in  equation  (83)  denotes  the  discrete  approximation  of  the  right-hand  side  of 
the  momentum  equations  (70)  at  the  node  (i,j,k): 

RHS  =  A,^-(-^  +  H  (84) 

34,  ■'  34, 

Also,  Atj.j.k  in  equation  (83)  is  the  time  increment  which,  for  reasons  discuss  later  on,  varies  in 
space  (local  time  stepping).  For  the  sake  of  convenience,  however,  in  the  rest  of  the  analysis  the 
(i,j,k)  subscript  has  been  dropped. 

VI.4  Pressure-velocity  coupling  algorithm 

The  system  of  discrete  governing  equations  (82)  and  (83)  can  not  be  integrated  in  time  in  its 
current  form  due  to  the  lack  of  an  evolution  equation  for  the  pressure  field.  The  discrete 
momentum  equation  (83),  however,  can  be  substituted  in  the  di.screte  continuity  equation  (82)  to 
obtain  a  Poisson  equation  for  the  pressure  field  at  the  intermediate  stage  ((-1).  As  discussed  in 
Sotiropoulos  (1991),  on  a  non-staggered  grid  layout,  the  so  resulting  pressure  equation  would  be 
exactly  equivalent  to  the  discrete  continuity  equation— since,  at  steady  state  the  pressure  equation 
reduces  to  the  continuity  equation  (82)— but  it  would  yield  oscillatory  solutions  for  the  pressure 
field  (odd-even  decoupling).  To  overcome  this  difficulty,  it  was  proposed  by  Sotiropoulos  (1991) 
to  derive  the  discrete  pressure  equation  starting  from  a  modified  form  of  the  discrete  continuity 
equation.  More  specifically,  the  proposed  discrete  "continuity  equation"  reads  as  follows: 
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(85) 


DIV  [(}(,,] 

where 

L[]  =  {  5^'(Jg“At5^'  )  +  5^-’(Jr2^t5^^  )  +  5^XJg^^At5^^  )}[] 

L[l={5.'(Jg‘’At5^'  )  +  5^XJr“At5.=  )  +  5^'(Jg’^At5^’  )}[] 

The  detailed  reasoning  behind  the  selection  of  the  source  term  in  the  right  hand  side  of  equation 
(85)  can  be  found  in  Sotiropoulos  (1991).  Herein,  it  suffices  to  note  that  the  .source  term  in 
equation  (85)— which  is  nece.ssary  to  guarantee  the  smoothness  of  the  computed  pressure  field-is 
proportional  to  the  difference  between  the  orthogonal  parts  of  two  discrete  approximations  of  the 
Laplace  pressure  operator;  the  one  that  results  by  discretizing  over  2A^  and  the  one  that  results  by 
discretizing  over  4A^.  The  positive  constant  y  is  introduced  to  control  the  size  of  the  source  term 
and  minimize  the  error  in  the  satisfaction  of  the  discrete  continuity  equation.  Numerical 
experimentation  with  a  variety  of  three-dimensional  flows  has  shown  that  values  of  y  between  0.01 
and  0. 1  are  sufficient  to  eliminate  the  odd-even  decoupling  of  the  pressure  nodes. 

By  incorporating  the  discrete  momentum  equation  (83)  into  the  discrete  "continuity 
equation"  (85),  the  following  equation  is  obtained  for  the  pressure  field  at  the  ((-1)  stage  (for  the 
sake  of  convenience  we  .set  m=^'- 1 ): 

( 1  -y)  L[P  "’I  +  y  L[P  "’]  +  N[P  "’]  =  -L  DIV[Q "]  -  ^  (8«) 

ai 

where 

N[P]  =  {5^’[JAt  (gi25^^+gi35^^)]  -H  5^XJAt  (g»25^'+g2354^)] 

+  5^XJAt  (g‘36^--Hg235^=)|}[PI 
and 

a  =  6F;(JAt^iX^’) 

where  ^  contains  discrete  approximations  of  the  convective  and  viscous  terms  which  appear  in  the 
^-momentum  equation. 


(89) 

(90) 


(86) 

(87) 
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In  Sotiropoulos  (1991),  the  pressure  equation  (88)  was  solved  using  the  point  successive 
relaxation  method.  In  the  present  work,  however,  in  order  to  accelerate  the  convergence  of  the 
pressure  equation  as  well  as  that  of  the  global  time  marching  procedure,  equation  (88)  is  solved 
using  the  Alternate  Direction  Implicit  (ADI)  method.  For  that  reason,  a  time  derivative  of  the 
pressure  is  introduced  in  equation  (88)  transforming  the  pressure-Poisson  equation  into  a 
diffusion-like  evolution  equation: 

-p  ^  +  (1-Y)  L[P'"1  +  Y  UPH  +  N[Pn  =  ■ .  •  (91) 

at 

where  p,  in  equation  (91),  is  a  positive  preconditioning  constant  introduced  to  accelerate  the 
convergence  to  steady  state  as  discussed  in  a  later  paragraph.  Incorporating  the  first-order  accurate 
Euler  implicit  temporal  linearization  scheme 

P  '^  =  P  m-i  +  ^  At  4-  0(  AO  =  P  "’-O  AP 
ot 

in  equation  (91),  one  obtains 

AP  -  4L  { ( 1  -Y)  L[AP]  -h  y  L[AP]  }  =  -  ^  PRHS  (92) 

|3  P 

where 

PRHS  =  (1-Y)  L[P  "’-^1  Y  L[P"’-']  +  NIP  -  -L  DIV[Q k  (93) 

ai 

Note  that  the  non-orthogonal  pressure  terms  (N[  ])  in  equation  (92)  are  being  treated  explicitly  in 
order  to  preserve  the  tridiagonal  character  of  the  resulting  system.  Also  the  dissipation  term  for  the 
pressure  equation  (source  term  in  equation  (85))  is  treated  implicitly  in  equation  (92).  Such 
treatment,  however,  would  result  (for  the  general  case  of  Y  1 )  to  a  pentadiagonal  system— in  the  i- 
sweep,  for  instance,  the  points  i-2,  i-l,i,  i+1  and  i4-2  would  be  involved.  To  avoid  the  inversion 
of  a  pentadiagonal  system,  the  L[  ]  operator,  which  introduces  the  i-2  and  i-t-2  off-diagonal  terms, 
is  eliminated  from  the  left-hand  side  of  equation  (92)  by  setting  the  y  parameter  equal  to  one-only, 
of  course,  in  the  left-hand  side  of  equation  (92).  Obviously,  this  step  has  no  effect  on  the  .steady 
state  solution  of  the  pressure  equation  since  its  left-hand  side  vanishes  at  convergence.  Application 
of  the  ADI  approximate  factorization  method  to  equation  (92)  gives: 
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[1  -^8p'(Jg‘iAt5p')][l  -^5,^(Jr2^tL)][l  AP 

P  s  <;  P  q  q  p  q  q 

=  - ^  PRHS 

P  (94) 


Equation  (94)  is  solved,  in  three  consecutive  sweeps,  using  the  Thomas  algorithm.  It  is  important 
to  point  out  that  the  implicit  solution  of  the  pressure  equation  (94)  increases  the  overall  CPU 
requirement  no  more  than  five  percent,  mainly  because:  i)  the  three  tridiagonal  matrices  to  be 
inverted  in  equation  (94)  depend  only  on  the  metrics  of  the  geometric  transformation  (the  three 
diagonals  of  each  matrix  can  be  computed  once,  at  the  beginning  of  the  calculation,  and  stored  for 
subsequent  use);  and  ii)  the  Thomas  algorithm  can  be  vectorized  by  solving  simultaneously  for  all 
the  points  on  a  plane  perpendicular  to  the  current  sweep  direction. 

It  is  well  known  that  the  ADI  factorization  is  unconditionally  stable  (Anderson  et  al.,  1984) 
when  applied  to  the  parabolic-in-time  heat  equation.  Of  course,  the  same  argument  can  not  be 
readily  extended  to  the  pressure  equation  without  rigorous  stability  analysis  of  the  system  of 
governing  equations.  Numerical  experimentation  has  shown,  however,  that  stable  solutions  for 
the  pressure  equation  (94)  can  be  obtained  for  time  steps  much  larger  than  those  used  for  the 
momentum  equations.  More  specifically,  values  of  the  preconditioning  constant  p  of  the  order  of 
0. 1  have  been  used  in  all  the  calculations  reported  in  the  following  sections.  This  allows  the 
pressure  equation  to  operate  with  an  effective  time  step  one  order  of  magnitude  higher  than  the 
momentum  equations,  and  this  results  in  significant  convergence  acceleration  of  the  overall  time 
marching  procedure. 

VI.5  Solution  of  the  k  and  £  equations 

The  k  and  e  equations  (74)  and  (75)  are  discretized  in  space  using  second-order  upwind 
differencing  for  the  convective  terms,  and  three-point  central  differencing  for  the  viscous  terms. 
Equations  (74)  and  (75)  are  of  the  same  type  as  the  momentum  equations  (parabolic  in  time,  elliptic 
in  space)  and,  thus,  the  Runge-Kutta  scheme  could  be  applied  to  integrate  them  in  time.  The  k-e 
equations,  however,  contain  source  terms  which  are  stiff-the  production  term  G,  for  instance,  is 
very  large  near  solid  boundaries  while  it  decays  rapidly  to  zero  in  the  outer  part  of  the  boundary 
layer.  The  use  of  an  explicit  time-marching  algorithm  to  integrate  a  stiff  set  of  equations  may 
impose  severe  time  step  limitations  and  have  overall  stability  problems.  To  avoid  such  undesirable 
complications,  equations  (74)  and  (75)  are  solved  by  employing  the  ADI  approximate  factorization 
scheme. 
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The  e-equation  (75)  is  used  to  demonstrate  the  linearizadon  and  approximate  factorization 
procedure  (a  similar  procedure  is  employed  for  the  k-equation).  Using  the  Euler  implicit  temporal 
discretization  scheme,  equation  (75)  is  linearized  in  time  as  follows: 

Ae  +  At  [  -Qi  ^  +  C£2  ~  5^‘  +  V-  8^‘ 

iC  K 


-  j  5^'  (JvV*  5^0  ]  Ae  =  -At  RHS^ 

where 


(95) 


RHS=V^d 


t'e  + 


V' 


5^‘e 


-4- 5.'  (JvVJ5p’8)-C^,fG+C,,|^ 


(96) 


(97) 


Re  CTg 


(98) 


The  source  terms  of  equation  (75)  have  been  linearized  in  time  as 

g::!iG  =  £^G  +  ^G 

k  k  k 

(99) 

(£!lllll  =  ^£!!}i  +  2e!LAp 

k  k  k 


In  order  to  preserve  the  tridiagonal  character  of  the  resulting  system,  the  non-orthogonal  terms  of 
the  diffusion  operator  have  been  treated  explicitly  in  equation  (95).  For  the  same  reason,  the 
convective  terms  in  the  left-hand  side  of  equation  (95)  have  been  discretized  using  first-order 
upwind  differencing.  This  less  accurate  spatial  discretization  does  not  have  any  effect  on  the 
steady  state  solution— the  left  hand  side  of  equation  (95)  vanishes  at  convergence— since  the 
standard  second-order  upwind  differencing  has  been  employed  in  the  right-hand  side.  Finally, 
applying  the  ADI  approximate  factorization,  equation  (95)  reads  as  follows: 


1  +  At 


1  +At 


V2+  6^2  +  yl-  .  i  5^2  (jv£g22  5^2jJ  j 

5^'  +  V^'  5^'  - 1 6^’  (jvV^  5^^)] j  Ae=  -At*  RHS 


(100) 


where 

* 

At 


(101) 


Equation  (100)  is  solved  using  a  vectorized  version  of  the  Thomas  algorithm,  in  a  similar  fashion 
as  described  in  the  previous  section  for  the  pressure  equation.  The  k-equation  is  linearized  and 
solved  in  an  identical  manner  as  described  above.  The  only  difference,  of  course,  being  that,  in  the 
context  of  the  two-layer  model,  the  k-equation  is  solved  all  the  way  to  the  wall  while  the  e-equation 
is  solved  only  in  the  outer  layer. 

Due  to  the  rapid  changes  in  the  cross-sectional  shape  and  area  of  some  of  the  duct 
configurations  to  be  considered  herein,  care  must  be  exercised  when  specifying  the  matching 
boundary  for  the  two-layer  k-e  turbulence  model.  In  all  the  previously  reported  calculations  with 
the  two-layer  model  (see,  for  instance,  Kim,  1991)  the  approach  followed  is  to  match  the  one-  with 
the  two-equation  model  along  a  pre-selected  coordinate  surface  which  is  broadly  located  within 
what  is  normally  the  logarithmic  region.  This  approach,  although  it  works  well  for  ducts  of 
constant  cross-section  (and  is  adopted  herein  for  similar  configurations),  is  obviously  not  adequate 
for  transition  ducts  where  the  cross-section  changes  from  circular  at  the  inlet  to  rectangular  at  the 
exit.  For  that  reason,  in  the  calculations  with  transition  ducts  we  choose  to  follow  a  more  rigorous 
approach,  that  is  to  match  the  two  models  at  points  where  the  turbulence  Reynolds  number  Rey  (= 
Vky/v)  is  approximately  equal  to  250  in  order  to  ensure  that  the  wall  damping  effects  are  negligible 
(Chen  and  Patel,  1988)  beyond  that  distance. 

VI.6  Summary  of  algorithm  and  convergence  acceleration  techniques 

Assuming  that  the  solution  at  the  "n"  time  level  is  known,  the  solution  at  the  "n-t-1"  time 
level  is  obtained  through  the  following  steps: 

1 .  Solve  the  k  and  e  equations  (74)  and  (75),  as  described  in  section  5,  and  calculate  the 
new  eddy-viscosity  field. 
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2.  Using  the  currently  known  pressure  field,  calculate  the  artificial  mass  source  term  for 
the  pressure  equation  (right-hand  side  of  equation  (85)).  To  save  computational  time 
the  artificial  dissipation  term  is  frozen  in  all  subsequent  Runge-Kutta  stages. 

3 .  For  (=1  to  4  (  Q^=  Q  for  f  =  4); 

(a)  Compute  the  right-hand  side  of  the  pressure  equation  as  given  by  equation  (93). 

(b)  Solve  the  pressure  equation  (94)  to  obtain  the  pressure  field  at  the  "f-1"  stage. 
Since,  the  steady  state  solution  is  of  interest,  only  one  ADI  iteration  is  performed  on 
the  pressure. 

(c)  Using  equation  (83)  compute  the  velocity  field  at  the  "/-stage"  and  return  to  step 
(a). 

4.  Repeat  steps  1  to  3,  until  convergence  is  reached. 

The  convergence  rate  of  the  time  marching  procedure  is  enhanced  by  employing  the  local 
time-stepping  technique  along  with  implicit  residual  smoothing.  The  time  increment  is  computed 
and  stored  for  every  node  as  follows  (see  Martinelli  (1987),  Kunz  and  Lakshminarayana  (1991)): 

At,  j  k  =  min  (At.  j  k^At"  j_k) 


where 


I 

At  =  CFL  min(  VI77,  /gTT,  ) 


(^  + V,)  (g^-Hg^^+g^^) 

In  the  above  equations,  CFL  and  Q.  denote  the  Courant-Friedrich-Lewis  number  (hyperbolic 
stability  criterion)  and  the  von  Neumann  number  (parabolic  stability  criterion),  respectively.  The 
CFL  number,  used  herein,  is  an  approximate  one,  since  it  is  based  only  on  the  local  length  .scales 
of  the  computational  grid.  Although  an  exact  CFL  number  should  involve  the  local  velocity  scales 
as  well,  we  chose  to  use  this  approximate  formulation  in  order  to  avoid  the  calculation  of  At^  at 
every  new  iteration  level— this  purely  geometric  variation  of  At*  has  been  found  adequate  on  highly 
stretched  meshes  (Sotiropoulos,  1991).  Typically,  the  parabolic  stability  con.straint  dominates  only 
in  the  near-wall  region  where  the  grid  spacing,  the  velocity  and  the  eddy  viscosity  approach  zero. 
Sufficiently  far  from  the  wall,  however,  the  hyperbolic  stability  criterion  dictates  the  choice  of  a 
stable  time  increment.  In  the  present  calculations,  the  selection  of  the  local  time  increment  based  on 
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both  the  hyperbolic  and  parabolic  stability  criteria  (equation  (102))  was  found  necessary  for 
stability  only  for  Reynolds  numbers  of  the  order  of  10^  or  greater.  At  lower  Reynolds  numbers, 
the  hyperbolic  stability  criterir  '  was  sufficient  for  obtaining  converged  solutions. 

The  implicit  residual  smoothing  was  first  proposed  by  Lerat  (see.  for  example,  Hollanders 
et  al.,  1985)  for  use  with  the  Lax-Wendroff  scheme  and  was  later  adopted  by  Jameson  (1983)  to 
accelerate  the  convergence  of  Runge-Kutta  .schemes.  In  the  present  study,  the  implicit  residual 
smoothing  is  applied  to  the  right-hand  side  of  the  momentum  equation  as  in  Sotiropoulos  (1991). 
More  specifically,  the  residual  calculated  in  equation  (84)  is  smoothed  by  the  constant  coefficient 
implicit  operator  to  define  a  new  residual: 

(1-  CD, &....)(!-  (D,a:.2)(l-  cdA3.3)  Ms'  =  RHS'’ 


where 


j.k 


^.j.k  ^  ^i-l.j.k 


A  T 

(A^  )“ 


The  constants  cOj,  co.,  and  co^  are  smoothing  parameters  which  are  of  the  order  of  one  and  their 

subscripts  indicate  that  they  can  be  chosen  differently  for  each  spatial  direction.  Equation  ( 104)  is 
solved  using  the  Thomas  algorithm  and  the  .smoothed  residual  replaces  the  residual  RHS  in 
equation  (83).  The  residual  smoothing  is  applied  at  every  stage  between  steps  (3.b)  and  (3.c)  in 
the  solution  procedure  described  in  the  beginning  of  this  section  .  The  implementation  of  the 
implicit  residual  .smoothing  in  the  four  stage  procedure  allows  the  use  of  higher  CFL  numbers  and 
consequently  leads  to  a  significant  acceleration  of  the  convergence  rate-Sotiropoulos  (1991) 
reported  a  fifty  percent  gain  in  convergence  speed.  Moreover,  the  vectorized  version  of  the 
Thomas  algorithm,  used  for  inverting  the  three  linear  operators  in  equation  (104),  increases  the 
overall  CPU  time  by  no  more  than  five  percent. 

The  smoothing  coefficients  in  equation  (104)  are  constant  in  each  spatial  direction  and, 
therefore,  one  can  expect  this  formulation  to  be  optimal  for  grids  that  are  not  highly  stretched. 
Surprisingly  though,  equation  (104)  has  worked  quite  well  (Sotiropoulos,  1991)  for  highly 
stretched  grids  with  large  aspect  ratios.  For  further  acceleration  of  the  convergence  rate,  however, 
Martinelli  (1987)  proposed  a  formulation— for  the  two-dimensional  compressible  Navier-Stokes 
equations— where  the  smoothing  coefficients  in  equation  (104)  are  functions  of  characteristic  wave 
speeds.  The  idea  behind  Martinelli's  suggestion  is  that,  since  the  minimum  local  grid  spacing 
dictates  the  maximum  allowable  local  time  step  for  stable  calculations,  more  smoothing  should  be 
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applied  in  the  direction  of  that  minimum  spacing.  If  the  same  smoothing  is  applied  in  the  other 
spatial  directions— where  the  grid  spacing  is  coarser  and  the  time  step,  as  computed  by  equation 
(104),  is  much  smaller  than  the  local  stability  limit— the  damping  properties  of  the  scheme  are 
impaired.  Martinelli’s  formulation  was  extended  to  three-dimensions  by  Radiespiel  et  al.  (1989) 
and  Liu  and  Jameson  (1992).  In  the  present  study  a  formulation  similar  to  that  of  Liu  and  Jameson 
is  adopted  as  follows: 
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In  the  above  equations,  CFL*  denotes  the  CFL  number  of  the  unsmoothed  scheme  while  CFL 
denotes  the  desirable  CFL  number.  It  can  be  easily  seen  that,  for  non-equal  grid  spacings,  larger 
smoothing  will  be  applied  in  the  direction  of  smaller  spacing,  while  for  uniform  grid  the  three 
smoothing  coefficients  are  equal.  In  addition,  these  equations  involve  only  geometric  quantities 
and,  consequently,  the  smoothing  coefficients  need  to  be  computed  only  once  at  the  beginning  of 
the  calculation. 


VII.  THE  DISCRETE  CONTINUITY  EQUATION  IN  METHODS  I  &  II 

It  is  well  known  that  on  a  non-staggered  computational  grid— which  is  used  by  both 
methods  employed  here— the  discrete  continuity  equation  cannot  be  satisfied  to  machine  zero  with 
the  resulting  pressure  field  being  smooth  (Strikwerda  (1984),  Sotiropoulos  (1991)).  In  general,  a 
smooth  pressure  field  can  be  obtained  only  at  the  expense  of  accuracy  of  the  discrete  continuity 
equation.  Thus,  a  successful  non-staggered  grid,  primitive  variable  method  must  have  built  into  it 
a  dissipative  mechanism  that  eliminates  the  pressure  decoupling  and,  at  the  same  time,  a 
mechanism  that  minimizes  any  errors  introduced  in  the  discrete  incompressibility  condition. 

Method  n  was  specifically  tailored  to  account  for  the  ambiguities  associated  with  the  non- 
staggered  grid.  The  artificial  mass  source  term,  introduced  in  the  discrete  continuity  equation  (85), 
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eliminates  the  decoupling  of  the  pressure  nodes  but  at  the  same  time  its  size  is  controlled  via  the 
constant  y.  An  estimate  for  the  size  of  this  source  term  can  be  easily  obtained  if  we  consider  its 
form  for  the  case  of  uniform  cartesian  grid  (see  Sotiropoulos,  1991): 

DIV[Q, .  J  =  -Yn^  +  •  •  ■  ]  (H)8) 


where  the  subscript  II  is  used  to  differentiate  from  an  equivalent  parameter  which  will  be 
subsequently  used  for  Method  1.  Also  the  dots  in  equation  (108)  imply  similar  terms  in  the  other 
spatial  directions,  which  are  omitted  for  convenience.  Recall  now  that  Method  II  employs  the  local 
time  stepping  technique  according  to  which  the  local  time  increment  is  proportional  to  the  grid 
spacing,  with  the  constant  of  proportionality  being  the  CFL  number  (At  =  CFL  Ax).  Incorporating 
this  in  equation  ( 108),  we  obtain 


3x4 


(109) 


This  shows  that  the  error  in  the  satisfaction  of  the  discrete  incompressibility  condition  is  third  order 
in  space  (lower  than  the  truncation  error  of  the  second-order  accurate  finite  difference 
approximations).  Also,  it  is  important  to  note  that  the  use  of  the  local  time  stepping  scales  the  time 
increment  out  of  equation  (109)  and  makes  the  steady  state  solution  independent  of  the  time  step. 
Moreover,  the  CFL  number  in  equation  (109)  does  not  have  any  significant  effect  in  the  steady 
state  solution,  since  its  maximum  allowable  value  for  stability  is  4  (CFL/4  <  1).  All  the  herein 
reported  calculations  performed  with  Method  II  use  Yu  =0.01^1.1. 

Let  us  examine  now  the  accuracy  of  Method  I  insofar  as  the  satisfaction  of  the  discrete 
continuity  is  concerned.  For  convenience,  but  without  loss  of  generality,  we  use  only  the 
derivatives  in  the  pressure  equation  (66).  First  recall  that  the  ^Cmomentum  equation  at  the  (i,j,k) 
node  reads  (equation  (64)): 


V,(1)  =  V,(1)-E‘‘® 

I, 

with  the  pressure  gradient  term  computed  as 


(110) 


iPj  -  P.-i 


(111) 


42 


The  pressure  equation  is  derived  from  the  discrete  continuity  equation 


by  substituting  in  this  equation  (110),  using  simple  linear  interpolation  to  compute  the 
pseudovelocities  at  the  cell  interfaces  (equation  (68)),  and  calculating  the  pressure  gradient  terms  as 
follows; 


h± 


P.±i  -  P. 


(112) 


The  resulting  pressure  equation  (66)  reads 


[Jg/{Mi)i...-[Jg7rv(i)],., 


_pii  P.^i-P.  pii 
“'^1+1/2  1  '^1-1/2  , 


(113) 


where  for  convenience  By  simply  adding  and  subtracting  appropriate  terms  and 

using  linear  interpolation  to  compute  the  coefficients  of  the  pressure  derivatives  at  the  interfaces, 
the  right-hand  side  of  the  pressure  equation  (113)  can  be  rearranged  as  follows: 

[Jg-,'(^V(l)l,„-[Jg,'pV(l)l,.,  ,  _  I  c'lP.-P-M 


-J—rH  ^  Pi-*-2~^Pn-l''~Pi  .  qI  lPi4^2~^P|'^Pi-l  ^  H  *  P|'^P|-l'''P|-2i  , 

‘+l  I  ‘  I  1-1  I  J 

2A^  2A^  2A^  2A^ 


or,  equivalently 


J-j- iag\'f),„[  V(1)W  -  E‘i, 


-(Jg-|'P)i.,[V(l)i.|-Ei',5j^||+. 
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(114) 


1  roll  ^i+2-2Pi+l  +  Pi  pllP,+2-2P,+  Pi.i  pll  P,-2P,.,  +  P,.2, 
fl'-i+l  1  1  +^,-1  1  1+ 

2A^  2A^  2A^  2A^ 

Note  that  equation  (114)  is  identical  to  the  original  pressure  equation  (113).  Assuming  that 
convergence  is  reached  for  the  pressure  and  the  momentum  equations,  equation  (110)  implies  that 
the  two  bold-faced  terms  in  the  left  side  of  equation  (114)  reduce  identically  to  the  velocity 
components  at  the  i+l  and  i-1  nodes.  Thus,  at  convergence,  the  pressure  equation  (114)  can  be 
written  as  follows: 

„1  2 

-V  KJg'i'fVd)).,,  -  (Jg;fV(l)),.,l+  -  =  55V  (C"55'5‘)P,  + . .  • 

2^^  (115) 
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The  right-hand  side  of  equation  (115)  is  the  artificial  mass  source  term  introduced  by  the  pressure-  ^ 
velocity  coupling  formulation  of  Method  I  in  the  continuity  equation.  This  artificial  source  term  is 
responsible  for  eliminating  the  odd-even  decoupling  of  the  pressure  nodes  on  the  non-staggered 
grid  arrangement.  It  is  interesting  to  point  out  that,  using  straightforward  algebraic  manipulations, 
the  artificial  source  term  (pressure  dissipation)  used  in  Method  II  (see  equation  (85))  can  be 
rearranged  (Sotiropoulos,  1991)  in  a  form  similar  to  the  right-hand  side  of  equat'on  (1 15)— the 
difference  being  that  the  Cll  coefficients  in  equation  (115)  depend  on  geometric  quantities  and  the 
velocity  field  (see  definition  of  EH)  while  the  corresponding  coefficients  in  Method  II  depend  on 
geometric  quantities  only.  Note,  however,  that  in  Method  II  the  source  term  (see  equation(85))  was 
added  explicitly  in  the  discrete  continuity  equation  before  deriving  the  equation  for  the  pressure.  In 
Method  I,  on  the  other  hand,  the  mass  source  is  implicitly  introduced  in  the  discrete  continuity 
equation  because  two  different  approaches  are  used  to  calculate  the  pseudovelocities  and  the 
pressure  gradient  terms  at  the  cell  interfaces.  Recall,  that  linear  interpolation  is  used  to  obtain  the 
velocities  at  the  cell  interfaces  (equations  (68))  while  the  pressure  gradient  is  discretized  using 
equation  (112).  If  the  same  interpolation  procedure  used  for  the  pseudovelocities  were  to  be  used 
for  the  pressure  gradient  terms,  no  error  would  be  introduced  in  the  discrete  continuity  equation  but 
the  resulting  pressure  equation  would  yield  oscillatory  solutions  for  the  pressure.  Assuming  now  a 
uniform  cartesian  grid,  equation  (115)  can  take  the  following  approximate  form: 
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where 


ei  = 


_ CpRe _ 

Ax  [1  +  Cp(C(j+C[))]+CpRe 


(117) 


Note  that  equation  (116)  is  only  an  approximation--as  opposed  to  equation  (108)  for  Method  II 
which  is  exact  on  a  uniform  cartesian  grid-because  the  part  of  the  C*  *  coefficient  that  contains  the 
finite  analytic  coefficients  (Ej  term)  has  been  locally  linearized  and  taken  out  of  the  derivative. 
Although  only  approximate,  equation  (i  16)  shows  that  Method  I  satisfies  the  discrete  equation  to 
accuracy  comparable  to  that  of  Methoo  II  (equation  (108)).  The  difference  between  the  two 
formulations,  however,  is  in  the  size  and  form  of  the  dissipation  coefficients  £j  and  Ej].  The 
coefficient  in  Method  I  (equation  (1 17))— and  consequently  the  computed  solution-is  a  function  of 
the  Reynolds  number  of  the  flow,  the  time  step  and  the  local  velocity  field  (a  dependency 
introduced  via  the  finite  analytic  coefficients  in  equation  (117)),  while  the  Eq  coefficient  in  equation 
(108)  is  a  user  specified  constant.  It  is  interesting  to  point  out,  however,  that  for  sufficiently  large 
Reynolds  numbers,  the  Ej  coefficient  is  of  the  order  of  one  (see  equation  (117)),  while  it 
approaches  zero  with  the  Reynolds  number.  It  can  be  deduced,  therefore,  that  in  the  general  case 
of  any  finite  Reynolds  number-and  assuming  that  the  left  term  in  the  denominator  of  equation 
(117)  is  always  positive--the  E[  term  is  positive  and  less  or  equal  to  one. 

VIII.  A  COMPUTATIONAL  COMPARISON  OF  METHODS  I  &  II 

In  this  section  a  computational  comparison  of  Methods  I  and  II  is  made  in  order  to  identify 
their  relative  merits  and  weaknesses,  and  decide  on  an  optimal  numerical  approach  best  suited  to 
the  problems  under  consideration.  For  this  purpose,  both  methods  are  employed  to  calculate  the 
following  three  cases  (see  also  table  1  for  more  details  on  each  case): 

(i)  laminar  flow  (Re  =  790)  through  a  90®  square  bend  with  fully-developed  entry  flow 
(Humphrey  et  al.,  1977), 

(ii)  laminar  flow  (Re  =  790)  through  a  90®  square  bend  with  developing’  entry  flow 
(Taylor  et  al.,  1982),  and 

(iii)  turbulent  flow  (Re  =  224,000)  through  a  90®  rectangular  duct  of  aspect  ratio  six  with 
developing  entry  flow  (Kim,  1992). 
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The  two  laminar  flow  cases— both  very  well  documented  in  experiments— are  selected  in  order  to 
assess  the  accuracy  of  the  numerics  without  the  uncertainties  of  a  turbulence  model.  The  third  test 
case,  on  the  other  hand,  is  selected  to  evaluate  the  combined  performance  of  the  numerics  with  the 
same  turbulence  model.  To  establish  a  basis  for  meaningful  and  fair  comparison,  „ii  sub^equendy 
reported  calculations  were  performed  with  identical  computational  meshes,  starting  from  the  same 
initial  conditions,  using  the  same  boundary  conditions,  and  with  the  same  convergence  criteria  tor 
both  methods.  The  relative  performance  of  the  methods  is  discussed  with  emphasis  on  numerical 
accuracy  and  overall  computational  efficiency. 

VIII.  1  Fully-developed  laminar  flow  through  a  90*^  square  bend 

The  measurements  of  Humphrey  et  ul.  (1977)  were  carried  out  at  Reynolds  number.  Re  = 
790.  with  a  corresponding  Dean  number,  De  =  368.  In  the  experiment,  a  long  straight  entry  duct 
was  used  to  realize  fully-developed  riow  at  the  entrance  of  the  bend  (see  table  1).  In  the  pre.sent 
calculations,  however,  the  solution  domain  starts  five  hydraulic  diameters  upstream  of  the  bend 
and  fully-developed  flow  conditions  are  prescribed  there  using  the  analytical  solution  given  in 
White  (1974). 

Calculations  are  carried  out  on  two  numerical  grids,  namely,  grid  A  with  74x41x21  nodes, 
and  grid  B  with  99x41x21  nodes,  in  the  streamwise,  radial  and  normal  directions,  respectively. 
The  streamwise  spacing  inside  the  bend  is  2°  for  grid  A  and  1.5°  for  grid  B.  The  grid  nodes  in  the 
cross-sectional  plane  are  distributed  using  a  hyperbolic  tangent  stretching  function  with  stretching 
ratios,  in  all  spatial  directions,  nowhere  exceeding  1 .3.  The  exit  boundary  for  both  grids  is  located 
seven  hydraulic  diameters  downstream  of  the  bend.  The  physical  and  computational  domains, 
along  with  the  coordinate  systems,  are  shown  in  figure  4.  Typical  views  of  the  grid  in  the 
symmetry  plane  and  the  cross-section  are  also  shown.  As  the  duct  geometry  is  symmetric  with 
respect  to  the  z-axis  and  as  the  entry  flow  profile  is  also  symmetric,  only  one-half  of  the  duct  was 
considered. 

As  mentioned  above,  the  same  set  of  boundary  conditions  is  employed  for  both  methods. 
More  specifically,  Dirichlet  conditions  are  used  for  the  velocity  components  at  the  inlet  (analytical 
solution)  and  on  the  solid  walls  (no-slip,  no-flux  condition),  while  at  the  exit  the  three  velocity 
components  are  computed  by  assuming  zero  streamwise  diffusion  (())^^  =  0).  The  pressure  at  the 
inlet,  exit  and  solid  boundaries  is  computed  using  linear  extrapolation  from  within  the  solution 
domain.  Finally,  on  the  symmetry  plane  (z  =  0),  the  governing  equations  are  solved  in  exactly  the 
same  way  as  for  every  internal  computational  node  using  mirror- image  reflection  for  the  grid  and 
the  flow  variables. 
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The  convergence  histories,  with  grid  A,  are  shown  in  figure  5  for  both  methods.  The 
vertical  axis  in  this  figure  is  the  logarithm  of  a  residual  defined  as; 
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where  N  is  the  total  number  of  grid  nodes,  A  denotes  changes  between  iterations  and  u,  v  and  w 
denote  the  physical  contravariant  velocity  components  for  Method  I  and  the  Cartesian  velocity 
components  for  Method  II.  The  horizontal  axis  in  figure  5  is  the  computational  work  expressed  in 
terms  of  Cray-2  CPU  minutes.  It  is  seen  that  Method  II  converges  twice  as  fast  as  Method  1.  This 
result,  however,  is  to  be  expected  because  Method  I  utilizes  the  full-transformation  approach  which 
increases  considerably  the  computational  work  per  grid  node  per  iteration.  The  increased  effort  is 
associated  with  discretization  of  the  additional  grid-related  terms  introduced  in  the  governing 
equations  to  account  for  the  spatial  variations  of  the  contravariant  base  vectors  (see  Section  V.  1). 
Recall  also  that  Method  0  is  a  mixed  explicit-implicit  method  and,  therefore,  the  computational 
work  required  per  grid  node  per  iteration  is,  by  default,  considerably  less  than  that  of  an  implicit 
method,  such  as  Method  I.  Moreover,  in  the  implicit  parts  of  Method  II  (implicit  residual 
smoothing  and  solution  of  the  pressure  equation)  the  operators  to  be  inverted  are  linear  and  thus 
they  are  computed  only  once,  at  the  beginning  of  the  iterative  procedure,  and  stored  for  subsequent 
use.  In  Method  I,  on  the  other  hand,  the  operators  to  t)e  inverted  (momentum  and  pressure 
equations)  are  nonlinear  and  they  need  to  be  re-computed  every  iteration. 

The  effect  of  grid  refinement  on  the  computed  solutions  is  shown  for  both  methods  in 
figi  re  6,  where  the  computed  streamwise  velocity  profiles— on  grid  A  (dashed  line)  and  grid  B 
(solid  line)— are  plotted  along  the  radial  direction  (from  the  inner  to  the  outer  wall)  for  several  axial 
locations  at  z  =  0.  and  z  =  -0.25.  The  comparison  between  grid  A  and  grid  B  solutions  is  very 
good  for  Method  I,  since  only  some  minor  changes  occur  at  0  =  60°  and  90°  near  the  duct 
centerline.  The  same  overall  trend  is  also  observed  in  the  solutions  obtained  by  Method  II,  except 
the  centerplane  profiles  at  0  =  60°  and  90°  where  a  local  maximum  of  the  velocity  appears  clearly 
near  the  inner  convex  wall  in  the  fine-grid  solution. 

The  fine-grid  solutions  of  the  two  methods  are  compared  with  the  experimental  data  of 
Humphrey  et  al.  (1977)  in  figures  7  and  8.  In  figure  7  the  computed  (dashed  line  for  Method  I  and 
solid  line  for  Method  II)  and  measured  (symbols)  streamwise  velocity  profiles  are  plotted  in  the 
same  fashion  as  in  figure  6.  The  calculations  are  in  good  agreement  with  each  other  and  with  the 
experimental  data  at  0  =  0°  and  30°.  Some  discrepancies  exist,  however,  at  the  more  downstream 
locations.  More  specifically,  at  0  =  60°  Method  I  yields  lower,  than  measured,  velocities  near  the 
bend  centerline,  while  at  0  =  90°  it  underpredicts  the  slope  of  the  streamwise  velocity  profile  at  the 
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z  =  -0.25  plane  near  the  inner  wall  and  fails  to  capture  the  local  velocity  maximum  near  the  inner 
wall  at  z  =  0.  The  computed  lower  velocities  near  the  bend  centerline  at  0  =  60°  are  even  more 
pronounced  in  the  solutions  obtained  with  Method  II.  Note  that  Method  II  predicts  two  well- 
defined  peaks  in  the  velocity  profile  at  z  =  -0.25  and  0,  in  contrast  with  the  Method  I  solutions 
where  two  velocity  peaks  appear  only  at  z  =  -0.25.  At  the  0  =  90°  plane,  on  the  other  hand. 
Method  II  correctly  predicts  the  slope  of  the  velocity  profile  near  the  inner  wall  at  z  =  -0.25  and 
captures  also  the  velocity  maximum  near  the  inner  wall  at  z  =  0. 

Figure  8  compares  the  computed  streamwise  velocity  profiles  and  the  corresponding 
experimental  data  at  0  =  30°,  60°  and  90°.  For  each  axial  location  the  velocity  profiles  are  plotted 
along  the  z-direction  (from  bottom  wall  to  symmetry  plane)  for  five  different  radial  locations.  As 
can  be  seen,  both  methods  yield  very  similar  solutions  at  0  =  30°  and  60°  but  some  discrepancies 
exist  between  experiment  and  calculations,  particularly  near  the  inner  wall.  At  0  =  90°  Method  II 
agrees  very  well  with  the  experimental  data,  while  Method  1  appears  to  grossly  underpredict  the 
streamwise  velocity  at  r*  =  0.9.  Figure  9,  which  is  of  the  .same  format  as  figure  8,  compares 
profiles  of  the  radial  velocity  component  calculated  by  the  two  methods;  no  experimental  data  for 
the  radial  velocity  component  are  available.  As  expected,  differences  between  the  two  calculations 
occur  at  the  same  locations  where  discrepancies  were  observed  in  the  streamwise  velocity 
component  (figure  8).  Finally,  a  more  global  picture  of  the  computed  solutions  can  be  obtained 

from  figure  10  where  contours  of  the  calculated  streamwise  velocity  component  are  compared  with 

* 

measurements  at  0  =  60°  and  90°. 

VIII.2  Developing  laminar  flow  through  a  90°  square  bend 

The  measurements  of  Taylor  et  al.  (1983)  were  made  at  Re  =  790  and  De  =  368.  The  duct 
and  bend  geometry  was  almost  identical  to  that  of  Humphrey  et  al.  (1977),  the  only  difference 
being  a  shorter  entry  length  in  order  to  obtain  a  thin  boundary  layer  at  the  entrance  to  the  bend. 
The  calculations  are  carried  out  on  a  single  computational  grid  with  69x41x21  nodes  in  the  axial, 
radial  and  normal  directions  (with  a  streamwise  spacing  inside  the  bend  of  3°).  No  grid  refinement 
study  is  carried  out  for  this  case  since  the  response  of  each  method  to  a  finer  grid  was  established 
in  the  previous  calculation.  It  should  be  pointed  out,  however,  that  the  grid  size  employed  here, 
particularly  the  number  of  planes  in  the  streamwise  direction,  is  probably  not  sufficient  for  grid 
independent  solutions  to  be  established. 

The  computations  started  7.5  hydraulic  diameters  upstream  of  the  bend  entrance  to  match 
the  experimental  configuration,  and  a  uniform  (plug)  inflow  velocity  profile  was  specified. 
Starting,  however,  with  uniform  flow  at  this  location  does  not  exactly  represent  the  experimental 
conditions  because,  in  reality,  a  boundary  layer  has  already  started  forming  at  that  location. 
Govindan  et  al.  (1991)  reported  that  in  order  to  match  the  experimentally  observed  thickness  of  the 
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boundary  layer  at  entrance  of  the  bend  they  had  to  start  their  calculations  8.5  diameters  upstream 
with  a  very  thin  boundary  layer.  No  such  effort  is  made  in  the  present  study  but,  as  the 
subsequently  presented  results  indicate,  the  effect  of  this  inaccurate  inlet  condition  is  small  and 
limited  to  the  first  one-third  of  the  bend  only.  The  exit  boundary  is  located  7  hydraulic  diameters 
downstream  of  the  bend.  The  boundary  conditions  are  applied  as  discussed  in  the  previous 
section. 

The  contours  of  streamwise  velocity  component  computed  by  the  two  methods  are 
compared  with  experimental  data  at  several  axial  locations,  inside  the  bend  and  in  the  upstream  and 
downstream  tangents,  in  figure  1 1.  The  effect  of  the  inlet  conditions  in  the  flow  development  can 
be  seen  primarily  at  x  =  -0.25dh  and  0  =  30°.  where  the  computed  velocities  in  the  region  near  the 

bend  centerline  are  lower  than  the  measured  ones.  This  trend-which  implies  that  the  core  flow  has 
not  accelerated  as  much  as  the  experimental  data  indicate~can  be  attributed  to  the  specification  of  a 
plug  flow  profile  at  the  inlet  of  the  computational  domain,  as  previously  discussed.  At  subsequent 
axial  locations,  the  distortion  of  the  measured  isovels  indicates  the  development  of  a  strong 
secondary  motion  which  increases  the  velocity  near  the  outer  wall  and  decreases  it  near  the  inner 
wall.  Method  II  appears  to  capture,  with  reasonable  accuracy,  the  major  features  of  the  flow  field. 
The  isovels  computed  by  Method  I.  on  the  other  hand,  are  not  distorted  to  the  extent  indicated  by 
the  experimental  data,  a  trend  which  implies  that  the  strength  of  the  secondary  motion  is  not 
correctly  predicted  (see,  for  example,  the  calculated  0.6  isovel  in  figure  1  Ic  which  intersects  the 
centerplane  at  the  same  point  as  the  measured  0.4  isovel). 

In  order  to  assess  the  results  of  the  calculations  in  more  detail,  the  computed  streamwise 
velocity  profiles  are  compared  with  the  experimental  data  at  0  =  30O,  60°,  77.5°  and  x  =  0.25dh 
and  2.5  d^,  in  figure  12.  For  every  axial  location  the  velocity  profiles  are  plotted  along  the  z- 

direction  (from  the  bottom  wall  to  the  symmetry  plane)  at  five  radial  locations.  Both  methods  yield 
satisfactory  results  at  0  =  30^^  and  60°.  At  the  next  two  downstream  locations,  however. 

significant  differences  appear  between  the  two  numerical  solutions.  More  .specifically.  Method  II 
agrees  very  well  with  the  experimental  data  except,  perhaps,  near  the  inner  wall  at  0.25dh  where 

the  velocity  is  underpredicted.  Method  I,  however,  significantly  underpredicts  the  velocity  near  the 
inner  wall  at  0  =  77.5°  and  at  x  =  ().25dh;  note  that  this  trend  was  also  present  in  the  fully- 

developed  entry  flow  case  at  approximately  the  same  locations.  Figure  13  compares  the  computed 
radial  velocity  profiles  with  the  experimental  data.  It  is  seen  that  the  two  numerical  solutions  are  in 
much  closer  agreement  in  this  case.  Differences,  however,  do  exist  at  the  same  locations  where  the 
disagreements  in  the  axial  velocity  component  al.so  occur.  For  instance.  Method  I  underpredicts 
the  crossflow  near  the  inner  wall  at  0  =  77.5°  while  it  overpredicts  it  at  0.25dh.  Also  Method  1 

fails  to  capture  the  experimentally  observed  structure  of  the  radial  velocity  profiles  (see  profiles  at 
r*  =  0.7  for  0  =  77.5°  and  for  0.25dh).  It  would  appear,  therefore,  that  the  observed  differences 
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in  the  predicted  strength  of  the  secondary  motion,  no  matter  how  small  they  may  be,  are 
responsible  for  the  significant  differences  in  the  distribution,  within  a  cross-section,  of  the 
streamwise  momentum. 

VIII.3  Turbulent  flow  in  a  90**  rectangular  bend  with  developing  entry  flow 

These  calculations  correspond  to  the  recent  experiment  of  Kim  (1992),  who  carried  out 
detailed  mean  flow  and  turbulence  measurements—as  well  as  numerical  calculations  using  Method 
I—  through  a  90®  rectangular  duct,  of  aspect  ratio  6,  at  Reynolds  number  Re  =  224,000  (based  on 
the  duct  width  and  the  mean  bulk  velocity).  An  overall  view  of  the  wind  tunnel  and  duct  geometry, 
as  well  as  the  sections  at  which  measurements  were  reported  (Ul,  U2,  15,  45,  75,  D1  and  D2), 
are  shown  in  figure  14.  As  can  be  seen  in  this  figure,  the  flow  enters  the  inlet  tangent  of  the 
curved  duct  through  a  short  transition  duct  (a  two-dimensional  6:1  contraction).  The  transverse 
pressure  gradients  on  the  top  wall  of  the  contraction  induce  a  pair  of  vortices  inside  the  top-wall 
boundary  layer  resulting  in  three-dimensional  flow  at  the  inlet  of  the  upstream  straight  tangent. 
Figure  15  shows  the  measured  mean  velocity,  turbulent  kinetic  energy  and  transverse  Reynolds 
stress  fields  at  station  Ul.  The  apparent  complexity  of  the  inlet  flow  requires  a  careful 
specification  of  inlet  conditions  for  the  numerical  calculations  in  order  to  properly  represent  the 
experimental  situation.  For  that  reason,  following  Kim  (1992),  the  experimental  data  at  station  U 1 
are  used  to  construct  appropriate  inlet  distributions  for  the  mean  velocity  components  and  the 
turbulent  quantities. 

The  computational  domain  starts  4.5H  upstream  from  the  inlet  of  the  bend  (section  U 1)  and 
extends  up  to  30H  downstream  from  the  exit  of  the  bend.  A  numerical  grid  with  62x69x52  nodes, 
in  the  streamwise,  radial  and  normal  directions,  respectively,  is  used  for  the  subsequently  reported 
calculations  with  both  Methods  I  and  II.  The  streamwise  spacing  inside  the  bend  is  while  the 
first  coordinate  surface  just  off  the  duct  walls  is  located  well  within  the  laminar  sublayer,  around 
0.75,  almost  everywhere.  Method  1  requires  approximately  2  hours  of  Cray-YMP  CPU  time 
to  reduce  the  residuals  by  three  orders  of  magnitude,  while  the  same  level  of  convergence  is 
achieved  in  75  minutes  of  CPU  time  by  Method  11. 

Measured  and  computed  profiles  of  the  three  mean  velocity  components  at  several  sections 
through  the  bend  and  in  the  upstream  and  downstream  tangents  are  shown  in  figure  16.  Both 
numerical  methods  yield  identical  results  at  section  15,  which  are  in  good  overall  agreement  with 
the  experimental  data.  At  the  next  downstream  station  (45)  the  two  computed  solutions  are  still  in 
very  close  agreement  with  each  other  but  discrepancies  between  experiment  and  calculations  appear 
near  the  outer  concave  wall  of  the  duct.  As  indicated  by  the  computed  velocity  profiles,  the 
predicted  boundary  layer  is  somewhat  thicker  than  the  measured  one.  This  discrepancy  should  be 
attributed  to  the  well  known  inability  of  k-e  based  models  to  capture  the  effects  of  concave 
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curvature— that  is  increased  levels  of  turbulent  mixing  which  result  in  fuller  velocity  profiles  near 
concave  walls-rather  than  to  numerical  inaccuracies  (Kim,  1992).  The  discrepancies  between 
experiment  and  calculations  in  the  prediction  of  the  boundary  layer  thickness  near  the  concave  wall 
persist  in  the  following  downstream  stations  75,  D1  and  D2  even  after  the  removal  of  the 
curvature.  The  agreement  between  the  two  computed  solutions  is  still  quite  close  at  sections  75 
and  D1  but  additional  discrepancies  between  experiment  and  calculations  appear  near  the  inner 
wall.  More  specifically,  as  indicated  by  the  shape  of  the  measured  velocity  profiles,  a  pressure- 
driven  longitudinal  vortex  develops  near  the  inner  wall  which  transports  low  momentum  fluid  from 
the  top  wall  boundary  layer  towards  the  inner  side  wall.  As  a  result,  the  measured  streamwise 
velocity  profiles  exhibit  an  S-shaped  pattern  as  seen  in  sections  75  (Z  =  0.5)  and  D1  (Z  =  0.5  and 
0.75).  Both  numerical  methods  fail  to  mimic  the  shape  of  the  measured  profiles,  a  trend  which 
implies  that  the  strength  of  the  secondary  motion  is  underpredicted.  This  can  be  seen  by  inspecting 
the  measured  and  computed  W-velocity  profiles.  Note,  however,  that  at  sections  75  and  D1 
Method  II  predicts  peak  values  of  the  W-velocity  (parallel  to  the  vertical  walls,  see  figure  14) 
component  which  are  consistently  higher— and  consequently  closer  to  the  measurements— than 
those  predicted  by  Method  I.  The  same  trend  in  the  prediction  of  the  secondary  motion  continues 
in  the  next  downstream  station  D2  which  is  the  last  station  at  which  measurements  were  taken.  At 
station  D2,  the  U-velocity  profiles  calculated  by  Method  II  exhibit  a  much  more  pronounced  S- 
shaped  structure  than  those  calculated  by  Method  I.  Naturally,  the  differences  in  the  computed 
streamwise  velocity  field  are  induced  by  the  differences  in  the  corresponding  transverse  velocity 
components,  since  Method  II  consistently  predicts  significantly  higher  peak  values  of  the  W- 
velocity  profile.  More  specifically,  at  Z  =  0.5  and  0.75,  Method  II  overpredicts  somewhat  the 
measured  W-component  but  good  agreement  between  measurement  and  calculation  is  achieved  at  Z 
=  1.00.  Method  I,  on  the  other  hand,  underpredicts  the  magnitude  of  the  secondary  motion  by  as 
much  as  50  percent  at  Z  =  1.0.  It  should  be  pointed  out,  however,  that  the  shape  and  magnitude  of 
the  U  and  W  velocity  profiles  at  Z  =  0.5  and  0.75  computed  by  Method  II  imply  that  the  predicted 
longitudinal  vortex  is  located  at  a  lower  position  (nearer  to  the  top  wall)  as  compared  to  the 
measured  one.  In  other  words.  Method  II  yields  a  better  prediction  for  the  strength  of  the 
secondary  motion  but  fails  to  accurately  capture  its  spatial  evolution. 

Figure  17  depicts  the  calculated,  by  .both  Methods  I  and  II,  and  measured  profiles  of  the 
turbulent  kinetic  energy  at  stations  45  and  Dl.  Both  methods  predict  very  similar  k-distributions 
which  are,  however,  in  gross  disagreement  with  the  measurements.  Particularly  near  the  concave 
wall,  the  calculations  consistently  underpredict  the  level  of  the  turbulent  kinetic  energy  by  as  much 
as  80  percent.  As  already  discussed  previously,  these  discrepancies  underline  the  inherent 
inadequacy  of  the  k-e  model  to  predict  the  high  levels  of  production  of  k  along  concave  walls. 
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Finally,  figures  18  and  19  show,  respectively,  the  computed  and  measured  contours  of 
streamwise  velocity  and  vorticity  at  several  locations.  These  figures  demonstrate  in  a  more  global 
sense  the  trends  discussed  in  the  previous  paragraphs  as  far  as  the  discrepancies  between 
experiments  and  calculations  as  well  as  those  between  the  two  numerical  methods  are  concerned. 
The  vorticity  contours  in  particular  (figure  19)  reveal  clearly  the  differences  between  the  two 
numerical  methods.  Although  both  methods  yield  similar  vorticity  fields  at  stations  75  and  Dl,  the 
vorticity  computed  by  Method  1  appears  to  diffuse  more  rapidly  in  the  next  downstream  station  D2. 
Note,  for  instance,  that  Method  Q  yields  a  core  vorticity  value  of  -0.9  which  is  more  than  twice  the 
corresponding  core  value  in  the  solution  obtained  by  Method  1. 

VIII.4  Discussion 

The  calculations  presented  in  the  previous  three  sections  reveal  significant  discrepancies 
between  the  steady-state  solutions  obtained  by  each  numerical  method,  despite  the  fact  that  identical 
computational  grids,  solution  domains  and  boundary  conditions  were  employed.  For  the  two 
laminar  flow  cases  computed,  it  is  seen  that  Method  1  tends  to  underpredict  the  magnitude  of  the 
streamwise  velocity  component  near  the  inner  (convex)  wall  for  bend  angles  larger  than  60*^  and 
that  this  trend  does  not  seem  to  be  affected  by  grid  refinement.  Method  11,  on  the  other  hand, 
appears  to  capture  correctly  most  of  the  flow  features  observed  in  the  experimental  data.  The 
differences  in  the  computed  axial  velocity  field  are  accompanied,  as  one  would  anticipate,  by 
qualitative  and  quantitative  differences  in  the  corresponding  computed  crossflow  velocities. 
Method  1,  for  instance,  tends  to  underpredict  the  strength  of  the  secondary  motion  inside  the  bend 
and  overpredict  it  near  the  end  of  the  bend  and  in  the  downstream  tangent. 

Similar  trends,  as  far  as  the  prediction  of  the  secondary  motion  is  concerned,  are  observed 
in  the  turbulent  flow  case  as  well  since  the  longitudinal  vortex  computed  by  Method  1  diffuses 
more  rapidly  in  the  downstream  straight  tangent.  However,  it  should  be  noted  that,  in  contrast 
with  the  laminar  flow  cases,  the  discrepancies  between  the  two  computed  solutions  become 
pronounced  only  several  duct  widths  downstream  from  the  end  of  the  bend.  Unfortunately,  it  is 
not  possible  at  this  stage  to  separate  the  effect  of  the  turbulence  model  from  that  of  the  numerical 
scheme.  Certain  deficiencies  of  the  two-layer  k-e  models— such  as  those  associated  with  the 
prediction  of  turbulence  along  concave  walls— have  already  been  discussed.  In  addition,  the 
present  study  (see  also  all  subsequently  reported  calculations)  confirms  the  results  of  Kim  (1992) 
and  W.  J.  Kim  (1992)  who  employed  the  same  turbulence  model  insofar  as  they  indicate  the 
inadequacy  of  isotropic  model  to  accurately  predict  the  origin,  growth  and  decay  of  complex 
vortical  flows.  Therefore,  it  is  concluded  that  the  turbulence  model  is  responsible  for  introducing  a 
large  amount  of  "false  eddy  viscosity"  into  the  numerical  solution,  and  this  apparently 
overshadows  the  effect  of  any  numerical  dissipation  and  discretization  errors  associated  with  the 
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numerical  methods.  This  conclusion  by  no  means  undermines  the  importance  of  an  accurate 
numerical  discretization.  Rather,  it  emphasizes  the  fact  that  the  numerics  become  crucially 
important  (as  our  laminar-flow  calculations  clearly  demonstrate)  when  the  turbulence  model  has 
been  sufficiently  refined  to  accurately  represent  real  flow  phenomena. 

The  discrepancies  between  the  two  numerical  solutions  are  certainly  more  dist  ■’’bing  than 
those  between  the  calculations  and  the  experimental  data—particularly  for  the  laminar  r'  .  ’  cases 
where  no  turbulence  model  related  uncertainties  are  present  in  the  calculations.  For  a  gi' 'tn  mesh 
size  (and  the  same  turbulence  model,  for  a  turbulent  flow  case)  one  would,  in  general,  anticipate 
differences  between  solutions  computed  by  two  numerical  methods,  mainly  due  to  differences  in 
the  order  of  the  truncation  error  inherent  in  the  numerical  discretization  of  continuous  differential 
operators.  Any  such  differences,  however,  should  decrease  with  increasing  grid  resolution  and 
eventually  approach  zero  in  the  continuum  limit,  provided  that  each  discretization  scheme  is 
consistent  (Peyret  and  Taylor,  1983).  Note  that  recent  calculations  carried  out  with  very  different 
numerical  methods  for  both  the  laminar  fully-developed  entry  flow  case  (Rogers  et  al.  ( 1991)  with 
an  artificial  compressibility  method  and  a  non-staggered  grid;  Rosenfeld  et  al.  (1991)  with  a 
fractional-step  method  and  staggered  mesh),  and  the  laminar  developing  entry  flow  case  (Williams, 
1991),  have  yielded  results  which  are  in  close  agreement  with  those  of  Method  II.  The  cause  for 
the  discrepancies  between  the  present  numerical  Methods  I  and  II  could  be  attributed  to  the 
different  approaches  each  method  adopts  to  formulate  and  discretize  the  flow  equations-such  as 
the  type  of  the  coordinate  transformation,  the  accuracy  in  the  satisfaction  of  the  discrete  continuity 
equation,  and  the  resolution  of  the  spatial  discretization  scheme.  These  factors  are  considered 
below. 

Generalized  coordinate  transformation  approach:  Method  1  utilizes  the  full-  transformation 
approach  with  the  physical  contravariant  components  as  independent  variables  as  opposed  to  the 
partial-transformation  approach  used  by  Method  II.  As  already  discussed,  the  full  transformation 
approach  introduces  additional  terms  in  the  governing  equations  to  account  for  the  spatial  variation 
of  the  coordinate  base  vectors  (Cristoffel  symbols)  and  imposes,  consequently,  additional 
smoothness  requirements  on  the  computational  grid.  As  a  result.  Method  I  is  expected  to  be  more 
sensitive  to  grid  discontinuities  which,  if  large,  could  deteriorate  the  accuracy  of  the  computed 
solution.  It  should  be  recalled,  however,  that  in  all  the  herein  reported  calculations:  i)  the 
configurations  under  consideration  are  geometrically  very  simple,  and  ii)  particular  care  has  been 
exercised  to  construct  smooth  numerical  grids  by  keeping  the  grid  stretching  ratio  less  than  1.3 
everywhere.  Moreover,  the  full-transformation  approach  was  also  adopted  by  Rosenfeld  et  al. 
(1991),  in  conjunction  with  a  finite-difference,  non-staggered  grid,  fractional-step  method,  but,  as 
mentioned  above,  their  results  (for  the  fully-developed  entry  flow  case)  are  in  very  good  agreement 
with  those  of  Method  II.  Therefore,  the  full-  transformation  approach,  although  it  may  have  some 
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impact  on  the  accuracy  of  Method  I,  would  not  appear  to  be  the  major  cause  for  the  observed 
differences  between  the  two  methods. 

Satisfaction  of  the  discrete  continuity  equation:  As  shown  in  section  V,  both  methods 
satisfy  the  discrete  continuity  equation  up  to  dissipative  terms  which  are  proportional  to  fourth 
order  derivatives  of  the  pressure.  There  are  two  major  differences  between  the  two  methods:  i)  the 
dissipative  term  in  Method  1  is  nonlinear  in  the  geometrical  quantities  as  well  as  the  velocity  field, 
while  in  Method  II  the  corresponding  term  is  nonlinear  only  in  the  geometrical  quantities;  and  ii) 
the  coefficient  of  the  dissipative  term  in  Method  I  (‘:ee  equations  (116)  and  (117))  is,  for 
sufficiently  high  Reynolds  number,  of  the  order  of  one,  while  the  corresponding  coefficient  used 
in  Method  II  is  a  user  specified  constant,  typically  set  equal  to  0.01.  These  differences  in  the 
accuracy  with  which  each  method  satisfies  the  discrete  continuity  equation  suggest  that,  for  a  given 
mesh  size,  differences  in  the  overall  accuracy  of  the  computed  solutions  are  to  be  expected. 
Sotiropoulos  (1991)  showed  that,  while  this  is  in  general  true,  accuracy  differences  between 
solutions  computed  with  different  values  of  the  Yjj  coefficient  in  equation  ( 108)  decrease  with  grid 

refinement.  This  is  also  to  be  expected  since  the  error  in  the  discrete  continuity  equation  for  both 
methods  is  proportional  to  the  grid  spacing.  In  the  present  calculations,  however,  a  34  percent 
increase  in  the  number  of  grid  nodes  (grid  B  for  the  first  laminar  flow  case  computed)  did  not  have 
any  significant  effect  on  the  solution  computed  by  Method  1  and  certainly  did  not  appear  to 
minimize  the  differences  between  the  two  methods. 

Spatial  discretization  scheme:  Method  I  utilizes  a  hybrid  finite-analytic  discretization 
scheme,  while  Method  II  uses  an  upwind  finite-difference  scheme  which  is  formally  second  order 
accurate  provided  the  grid  stretching  ratios  are  kept  close  to  unity  everywhere.  The  accuracy  of  the 
hybrid  finite-analytic  scheme,  on  the  other  hand,  can  not  be  readily  e.stimated— using  standard 
techniques  such  as  Taylor  series  expansion— due  to  the  complexity  of  the  finite-analytic 
discretization  coefficients  which  involve  hyperbolic  functions  and  exponentials.  Note,  however, 
that  the  hybrid  finite-analytic  scheme  employed  here  is  a  simplified  (for  reasons  of  computational 
efficiency)  9-point  version  of  the  more  general— and  consequently  more  ac''urate-28-point  scheme. 
This  9-point  scheme  has  worked  quite  well  for  three-dimensional  flows  where  the  streamwise 
direction  is  the  dominant  flow  direction  (Patel  et  al.,  1988;  Kim,  1991)  but  it  may  be  inadequate  for 
highly  swirling  flows  such  as  those  considered  in  the  present  work.  In  the  laminar  fully- 
developed  entry  flow  case,  for  instance,  there  are  regions  of  the  flowfield  where  the  radial  velocity 
component  is  as  high  as  50-60  percent  of  the  streamwise  component.  In  a  very  recent  work,  Yeo 
et  al.  (1991)  calculated  the  duct  flow  of  Humphrey  et  al.  (1977),  the  first  case  computed  in  the 
present  study  (section  VIII.  1),  using  three  different  finite-difference  schemes  for  discretization  of 
the  convective  terms,  namely,  the  first-order  upwind  differencing,  the  second-order  upwind 
differencing,  and  the  QUICK  scheme.  Interestingly  enough,  a  comparison  of  their  first-order 
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upwind  solutions  with  the  present  Method  I  (finite-analytic  discretization)  solutions  reveals  very 
similar  trends  in  the  region  between  the  middle  of  the  cross-section  and  the  inner  wall,  although 
Method  I  is  more  accurate  near  the  outer  wall.  The  similarities  and  differences  between  the  first- 
order  upwind  solution  of  Yeo  et  al.  and  the  present  Method  1  solution  would  suggest  that  the 
hybrid  finite-analytic  discretization  switches  to  first-order  upwinding  near  the  inner  wall  while  it 
maintains  higher  accuracy  elsewhere.  The  switch  to  first-order  upwind  could  be  probably 
attributed  to  high  values  of  the  local  cell  Reynolds  number  caused  by  the  very  large  secondary 
velocity  components— recall  that  the  finite-analytic  method  automatically  becomes  one  sided  for 
large  cell  Reynolds  numbers.  Another  interesting  point  regarding  the  calculations  of  Yeo  et  al.— 
which  is  not  addressed  at  all  by  the  authors-is  the  fact  that,  as  the  grid  is  refined,  the  second-order 
upwind  and  QUICK  solutions  tend  to  approach  each  other,  as  it  should  be  the  case  for  consistent 
finite-difference  schemes,  while  the  first-order  upwind  solution  improves  only  near  the  outer  wall 
and  maintains  essentially  the  same  trends  near  the  inner  wall.  Similarly,  the  grid  refinement  in  the 
present  calculations  did  not  seem  to  bring  Methods  I  and  11  closer  to  each  other  (see  figure  6).  An 
explanation  for  this  rather  odd  and  unexpected  behavior  may  be  found  in  the  recent  work  of  Brandt 
and  Yaneh  (1991).  They  proved,  mathematically  and  computationally,  that  when  the  grid  lines  are 
not  aligned  with  the  streamlines,  which  is  usually  the  case  except  for  some  very  cases  such  as 
boundary-layer  flows,  the  first-order  upwind  Difference  is  inadequate  regardless  of  the  grid 
spacing.  To  be  more  precise,  their  analysis  assumes  that  the  grid  spacing  tends  to  zero  but  it 
always  remains  bounded  from  below  by  a  certain  relation  which  involves  the  Reynolds  number 
(cell  Reynolds  number  constraint)  so  that  the  upwind  differencing  is  necessary  for  stability;  for 
sufficiently  small  grid  spacing  no  such  problems  arise  since  centered  schemes  can  be  employed 
without  stability  problems.  The  use  of  the  first-order  upwind  scheme  to  discretize  the  convective 
terms  introduces  nonisotropic  artificial  viscous  terms  which,  even  when  their  coefficients  tend  to 
zero,  may  effect  the  solution  significantly. 

In  summary,  it  appears  that  the  discrepancies  between  the  two  numerical  methods  used  here 
may  be  primarily  attributed  to  the  different  discretization  schemes  that  are  employed.  Improvements 
in  the  spatial  resolution  of  Method  I  could  be  achieved  by  replacing  the  9-point  hybrid  finite- 
analytic  scheme  with  a  more  general  one  involving  perhaps  the  corner  points  at  the  upstream  and 
downstream  faces  of  the  solution  element. 
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IX.  SOME  ADDITIONAL  CALCULATIONS  OF  FLOW  IN  DUCTS  OF 

SIMPLE  CROSS-SECTIONS 


In  this  section  we  report  a  series  of  laminar  and  turbulent  flow  calculations  for  flow  in 
curved  ducts  of  square  and  circular  cross-sections.  All  calculations  were  performed  with  Method 
n.  The  following  cases  are  computed: 

i)  laminar  flow  through  a  strongly  curved  90*^  pipe  bend,  with  fully-developed  entry  flow, 
measured  by  Bovendeerd  et  al.  (1987); 

ii)  laminar  flow  through  a  strongly  curved  90°  pipe  bend,  with  developing  entry  flow, 
measured  by  Enayet  et  al.  (1983); 

iii)  turbulent  flow  through  a  mildly  curved  180°  pipe  bend,  with  fully-developed  entry  flow, 
measured  by  Rowe  (1970); 

iv )  turbulent  flow  through  a  strongly  curved  1 80°  pipe  bend,  with  fully-developed  entry 
flow,  measured  by  Azzola  and  Humphrey  (1984); 

v)  turbulent  flow  through  a  strongly  curved  90*^  square  bend,  with  fully-developed  entry 
flow,  measured  by  Humphrey  et  al.  (1981);  and 

vi)  turbulent  flow  through  a  strongly  curved  180°  square  bend,  with  fully-developed  entry 

flow,  measured  by  Chang  (1983).  ‘ 

IX.l  90®  pipe  bend  of  Bovendeerd  et  al.  (1987) 

The  measurements  of  Bovendeerd  et  al.  (1987)  were  carried  out  at  Reynolds  number.  Re  = 
700  (based  on  the  bulk  velocity  and  the  diameter  D  of  the  pipe)  and  Dean  number,  De  =  286.  The 
radius  of  curvature  of  the  bend  is  3D,  and  the  length  of  the  upstream  tangent  is  50D  to  ensure  fully- 
developed  flow  at  the  entrance  of  the  bend. 

In  the  present  study,  the  calculations  start  5D  upstream  of  the  bend  where  a  fully- 
developed  velocity  profile  is  specified.  The  exit  boundary  is  located  7D  downstream  the  end  of  the 
bend.  The  grid  topology  and  the  related  nomenclature  are  depicted  in  figure  20.  Note  that  the  bend 
is  symmetric  with  respect  to  the  z-axis  and  since  the  inlet  flow  is  also  sym  ric,  only  one-half  of 
the  bend  needs  to  be  computed.  The  curvilinear  coordinates  (^,ti,Q  are  aligned  with  the  centerline, 
the  circumferential  direction,  and  the  tangential  direction,  respectively.  This  choice  of  coordinates 
facilitates  the  description  of  the  pipe  bend  geometry  at  the  expense  of  making  all  the  points  on  the 
bend  axis  singular.  In  the  transformed  domain,  the  bend  centerline  is  mapped  into  a  plane 
(ABCD),  as  shown  in  figure  20.  The  boundary  conditions  at  the  exit  (GLCB  plane),  the  solid  wall 
(FGKE  plane)  and  the  symmetry  boundaries  (AFGB  and  DEKC  planes)  are  applied  as  described 
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for  the  previous  square-bend  calculations,  while  a  simple  averaging  is  used  to  obtain  all  dependent 
variables  on  the  singular  plane  ABCD.  Note  that  the  grid  topology  depicted  in  figure  20  is  used 
through  out  in  all  calculations  with  pipe  bends.  The  results  for  this  case  are  obtained  with 
69x33x31  grid  nodes  in  the  axial,  circumferential  and  radial  directions;  finer  grid  calculations 
yielded  solutions  almost  identical  with  this  grid.  Typical  views  of  the  algebraically  generated  grid 
are  shown  in  figure  21. 

The  computed  contours  of  the  streamwise  velocity  component  are  compared  with  the 
experimental  data,  at  several  planes  inside  the  bend,  in  figure  22.  It  is  seen  that  the  calculations 
reproduce  with  remarkable  accuracy  the  development  of  the  streamwise  flow  as  it  passes  through 
the  bend.  In  figure  23,  the  computed  streamwise  velocity  profiles  on  the  symmetry  plane  are 
compared  with  the  experimental  data  at  several  streamwise  locations.  Again,  excellent  agreement  is 
observed  between  calculation  and  experiment. 

IX.2  90“  pipe  bend  of  Enayet  et  al.  (1983) 

Enayet  et  al.  (1983)  carried  out  measurements  at  Re  =  1096  and  De  =  693.  The  radius  of 
curvature  of  their  bend  was  2.8D  and  the  length  of  the  upstream  tangent  was  only  5D  so  that  the 
flow  at  the  entrance  of  the  bend  was  still  developing. 

The  calculations  start  5D  upstream  of  the  entrance  to  the  bend  with  a  uniform  inlet  velocity 
distribution.  The  exit  boundary  is  located  7D  downstream  of  the  exit  from  the  bend.  Grid 
refinement  studies  showed  that  a  grid  with  89x43x41  points  in  the  streamwise,  circumferential  and 
radial  directions,  respectively,  is  necessary  to  obtain  grid  independent  solutions.  Note  that,  for  the 
pipe  bend  of  Bovendeerd  et  al.,  a  much  coarser  grid  (69x33x31)  yielded  a  grid  independent 
solution.  This  is  attributed  to  the  fact  that,  in  the  present  case,  the  boundary  layer  entering  the  bend 
is  much  thinner— as  compared  to  the  fully-  developed  profile  in  the  previous  case— and 
consequently,  finer  grids  are  required  to  accurately  resolve  the  steep  velocity  gradients  in  the  near¬ 
wall  region. 

The  calculated  contours  of  streamwise  velocity  are  compared  with  the  experimental  data  at 
four  streamwise  locations  in  figure  24.  At  the  9  =  30“,  60“  and  75“  cross-sections,  the 
calculations  reproduce  quite  well  the  overall  features  of  the  flowfield.  Some  discrepancies 
observed  between  experiment  and  calculations,  primarily  in  the  high  velocity  region,  are  probably 
caused  by  the  specification  of  the  plug-flow  inlet  profile  at  a  location  where  the  boundary  layer  in 
the  experiment  has  already  started  to  develop.  As  the  secondary  motion  takes  over  the  dynamics 
of  the  flow  development  inside  the  bend,  however,  the  effect  of  the  inlet  flow  profile  is 
significantly  reduced.  This  can  be  seen  from  the  results  at  the  last  downstream  station  where  the 
calculations  reproduce  the  experimental  data  with  remarkable  accuracy. 
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IX.3  180®  pipe  bend  of  Rowe  (1970) 

Rowe  (1970)  carried  out  measurements  in  turbulent  flow,  at  Re=  236,000  and  De  = 
68,127,  in  a  180®  mildly-curved  pipe  bend,  with  radius  of  curvature  of  12D.  The  upstream 
tangent  was  69D  long  and  sufficiently  rough  to  ensure  that  the  pipe  inlet  had  an  axisymmetric  and 
fully-developed  turbulent  velocity  profile. 

A  grid  with  96x33x41  points  in  the  axial,  circumferential  and  radial  directions, 
respectively,  is  used  for  the  present  calculations.  Grid  dependence  studies  with  coarser  grids 
indicate  that  this  grid  is  sufficient  for  grid  independent  solutions.  The  first  =  constant  coordinate 
surface  of  the  solid  wall  is  located  at  y"^  =  1  almost  everywhere.  To  obtain  the  necessary  inlet 
conditions,  the  flow  through  a  very  long  straight  pipe  is  calculated  and  the  resulting  fully- 
developed  profiles  of  velocity  and  turbulence  parameters  are  specified  one  pipe  diameter  upstream 
of  the  entrance  to  the  bend. 

Figure  25  depicts  the  computed  and  measured  contours  of  the  velocity  head,  at 

1  T 

several  cross-sections  within  the  bend— the  velocity  head  is  nondimensionalized  with 

where  Uq  is  the  centerline  velocity  at  the  inlet.  At  the  0  =  0®  plane,  the  measured  velocity 
contours— which  should  correspond  to  fully  developed  turbulent  pipe  flow— imply  smaller  than  the 
calculated  velocities  in  the  core  region.  This  discrepancy  indicates  that  the  flow  in  the  experiment 
was  still  developing  at  the  entrance  of  the  bend  (such  a  conclusion  is  also  supported  by  the 
calculations  presented  in  the  next  section,  where  a  straight  pipe  calculation  is  also  carried  out, 
although  at  a  lower  Reynolds  number,  and  the  computed  fully-developed  axisymmetric  velocity 
profile  is  in  very  good  agreement  with  the  measurements).  The  inlet  velocity  profile  affects  the 
predicted  streamwise  flow  development  inside  the  bend,  since  the  velocities  in  the  core  are 
consistently  higher  than  the  measurements.  Overall,  however,  the  calculations  predict  fairly  well 
the  shifting  of  the  maximum  of  the  velocity  towards  the  outer  radii  and  the  accumulation  of  low 
speed  fluid  near  the  inner  radii.  But  the  latter  process  is  not  as  pronounced  in  the  calculations,  as 
the  experimental  data  indicates,  since  higher  velocities  are  predicted  near  the  inner  radii.  Although 
no  crossflow  measurements  were  reported  by  Rowe,  we  can  speculate  that  this  trend  is  due  to  an 
underestimation  of  the  strength  of  the  secondary  motion  which  is  re.sponsible  for  the  distribution  of 
the  streamwise  momentum  within  the  cross-section. 

IX.4  180®  pipe  bend  of  Azzola  and  Humphrey  (1984) 

Azzola  and  Humphrey  measured  the  flow  through  a  1 80®  strongly  curved  pipe  bend  (R^;  = 

3.375D)  at  two  Reynolds  numbers.  Re  =  57,400  and  110,000,  corresponding  to  De  =  31,300  and 
59,900.  The  entrance  length  of  their  bend  was  54.7D  to  ensure  axisymmetric  fully-developed  flow 
at  the  inlet  to  the  bend. 
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Calculations  are  carried  out  only  for  the  low  Reynolds  number  case.  The  computational 
grid  comprises  154x41x23  nodes  in  the  streamwise,  radial  and  circum-ferential  directions, 
respectively.  The  first  ^  =  constant  coordinate  surface  just  off  the  pipe  wall  is  located  at  y+  =  0.6 
almost  everywhere.  The  inlet  profiles,  for  the  velocity  and  turbulent  quantities,  are  obtained  by 
solving  the  developing  flow  through  a  long  straight  pipe.  The  calculations  start  two  diameters 
upstream  of  the  entrance  to  the  bend  and  continued  until  eight  diameters  downstream. 

The  computed  profiles  of  the  longitudinal  (ue)  and  circumferential  (U(i))  velocity 
components  are  compared  with  the  measurements  at  several  streamwise  locations  in  figure  26  (the 
profiles  are  plotted  from  the  symmetry  plane  to  the  wall  along  the  radius  located  90o  from  the 
symmetry  plane).  The  calculations  capture  well  the  development  of  the  longitudinal  velocity 
component  inside  the  bend  and  in  the  downstream  tangent— some  differences  exist  only  at  0  =  90^ 
where  the  measurements  indicate  a  larger  dip  in  the  velocity  profile  near  the  symmetry  plane.  This 
is  not  the  case,  however,  with  the  profiles  of  the  circumferential  velocity  component.  The 
calculations  underpredict  the  secondary  motion  in  the  first  quarter  of  the  bend  with  the  maximum 
discrepancy  occurring  at  45®  in  the  near-wall  region.  Further  downstream,  the  measurements 
indicate  a  reversal  of  the  secondary  motion  near  the  symmetry  plane  at  approximately  90“.  This 
process  is  predicted  qualitatively  by  the  calculations,  but  it  is  weaker  and  somewhat  delayed  since  it 
occurs  after  90®.  Moreover,  the  calculations  indicate  two  more  crossflow  reversals  occurring  at 
177^  and  at  one  diameter  downstream  of  the  end  of  the  bend.  The  measured  crossflow,  on  the 
other  hand,  exhibits  a  trend  to  reverse  direction,  decreases  between  135°  and  177®  and  increases 
between  177^  and  one  diameter  downstream  the  bend,  but  remains  always  positive  near  the 
symmetry  plane.  The  two  additional  crossflow  reversals  near  the  symmetry  plane  were  also 
predicted  (at  the  same  locations  as  the  present  calculations)  by  Azzola  et  al.  (1986),  who  calculated 
the  same  bend  on  a  mesh  of  similar  size  using  the  standard  k-e  model. 

The  effect  of  the  crossflow  on  the  streamwise  flow  development  can  be  seen  in  figures  27, 
where  the  calculated  contours  of  streamwise  velocity  component  and  the  cross-flow  velocity 
vectors  are  plotted  at  several  longitudinal  locations  inside  the  bend  and  in  the  downstream  tangent. 
At  90°  and  135°,  for  instance,  the  strong  secondary  motion  has  displaced  the  maximum  of  the 
velocity  in  the  interior  of  the  cross-section,  while  the  effect  of  the  crossflow  reversal  near  the 
symmetry  plane  is  clearly  visible  at  177°  and  x=-t-l  (see,  for  example,  the  change  in  the  shape  of 
the  U  =  1.0  isovel  between  these  two  locations). 

IX.5  90°  square  bend  of  Humphrey  et  al.  (1981) 

Humphrey  et  al.  (1981)  carried  out  turbulent  flow  measurements  in  a  square  bend 
geometrically  identical  to  that  used  by  Humphrey  et  al.  (1977)  and  Taylor  et  al.  (1982).  Their 
measurements  were  made  at  Re  =  40,000,  corresponding  to  De  =  26,000.  The  length  of  the  entry 
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tangent  was  45  hydraulic  diameters  so  that  fully-developed  flow  was  expected  at  the  entrance  of  the 
bend. 

This  experiment  of  Humphrey  et  al.  (198’ )  was  one  of  the  test  cases  selected  for  the  1980- 
81  Stanford  Conference  on  Complex  Turbulent  Flows  (Case  0512)  (Kline  et  al.,  1981)  and  for 
which  numerical  solutions  were  submitted  by  six  computer  groups.  The  turbulence  models  used 
for  the  calculations  by  the  various  computers  included:  i)  an  algebraic  eddy-viscosity  model  (one 
computer);  ii)  the  standard  k-e  model  with  wall  functions  (three  computers);  iii)  a  k-1  model  with 
wall  functions  (one  computer);  and  iv)  an  algebraic  Reynolds-stress  model  (one  computer).  These 
calculations  were  started  7.5  hydraulic  diameters  upstream  of  the  entrance  of  the  bend  where  the 
experimental  data  of  Melling  and  Whitelaw  were  used  to  obtain  inlet  profiles  for  the  flow  variables. 
A  sample  of  the  calculations  is  shown  in  figure  28  (taken  from  Vol.  II,  pp.  942-943,  of  the 
Conference  Proceedings).  Figure  28a  shows  general  agreement  on  the  streamwise  velocity 
component  at  the  inlet  plane  but  significant  disagreement  at  the  exit.  On  the  other  hand,  the 
predictions  of  the  radial  velocity  component  (Figure  28b)  are,  with  one  exception,  reasonable  at 
both  the  inlet  and  exit  planes.  Overall,  however,  none  of  the  six  computers  produced  results  that 
could  be  characterized  as  satisfactory.  Note  that  in  assessing  the  performance  of  the  numerical 
methods  applied  to  incompressible  duct  flows  in  the  Stanford  Conference,  J.  B.  Jones  concluded 
that  (Vol.  11,  pp.  914-918);  "among  all  the  cases  involving  secondary  flow  of  the  second  kind. 

Case  0512  has  produced  the  least  satisfactory  results." 

* 

The  present  calculations  are  carried  out  on  a  54x49x26  mesh  which  is  found  to  be  adequate 
for  grid  independent  solutions.  The  inlet  conditions  are  obtained  by  performing  a  developing  flow 
calculation  through  a  long,  straight,  square  duct.  The  fully-developed  solution  is  then  used  to  set 
the  inlet  profiles  of  the  velocity  components  and  the  turbulent  quantities,  2.5  hydraulic  diameters 
upstream  of  the  inlet  of  the  bend;  the  experimental  data  of  Humphrey  et  al.  (1983)  indicate  that  the 
flow  at  that  location  corresponds  to  fully-developed  square  duct  flow  with  no  influence  of  the 
downstream  bend.  The  calculated  fully-developed  mean  velocity  profile  is  compared  with  the 
measurements  of  Humphrey  et  al.  (at  x=-2.5  dh)  in  figure  29.  The  calculated  isovels  do  not  mimic 
the  shape  of  the  measured  ones  which  appear  to  be  bulging  towards  the  corners  of  the  cross- 
section.  Recall,  however,  that  the  distortion  of  the  measured  isovels  is  caused  by  the  existence  of 
secondary  motion  of  the  second  kind  which  is  induced  by  the  anisotropy  of  the  Reynolds  stresses. 
The  inability  of  the  calculations  to  reproduce  this  behavior  is  due  to  the  isotropic  eddy-viscosity 
model  employed  in  the  present  calculations.  The  present  model  implies  equal  normal  Reynolds 
stresses  and,  therefore,  cannot  predict  the  stress-driven  .secondary  motion. 

The  computed  and  measured  contours  of  streamwise  mean-velocity  component,  at  several 
cross-sections  within  the  bend,  are  shown  in  figure  30.  At  the  inlet  of  the  bend,  the  calculations 
indicate  that  the  high  speed  core  of  the  flow  shifts  towards  the  inner  wall.  This  trend  is  broadly 
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consistent  with  the  data  but  the  measured  velocities  in  the  high  speed  core  are  higher  than  the 
calculated  and  also  the  velocity  maximum  stays  closer  to  the  centerline  of  the  bend.  Moreover,  as 
indicated  by  the  shape  of  the  measured  isovels  near  the  corners,  the  effect  of  the  stress-driven 
secondary  motion  is  still  present  in  the  measurements  at  that  location.  At  the  0  =  45<’  cross- 
section,  the  calculations  imply  that  the  high  speed  core  continues  to  move  towards  the  inner  wall 
while  in  the  measurements  this  trend  has  been  reversed,  since  the  velocity  maximum  appears  to  be 
shifting  towards  the  centerline.  At  the  two  subsequent  cross-sections  the  calculations  predict  in 
general,  although  not  to  the  extent  indicated  by  the  data,  the  accumulation  of  low  speed  fluid  near 
the  inner  wall  and  the  shifting  of  the  high  speed  core  towards  the  outer  wall.  But  the  predicted 
velocities  in  the  high  speed  core  are  significantly  lower  than  the  measured.  Finally,  the  predicted 
and  measured  contours  of  radial  velocity  component  at  0  =  9(P  are  shown  in  figure  31.  The 
computed  radial  velocity  contours  are  in  good  agreement  with  the  measurements  except  near  the 
endwall  of  the  duct,  where  the  calculations  predict  .somewhat  higher  crossllow  velocities. 

The  reason  for  the  significant  disagreement  between  the  calculations  and  measurements  is 
not  altogether  clear  although,  as  noted  earlier  in  the  case  of  the  rectangular  duct  flow  (section 
VIII.3),  the  turbulence  model  mu.st  be  suspected.  To  some  extent,  and  particularly  near  the  inlet  of 
the  bend,  the  disagreement  can  be  attributed  to  the  simple  inlet  conditions  used  in  the  present 
calculations.  The  flow  through  the  bend  undergoes  very  rapid  changes  over  a  short  distance  (the 
bend  under  consideration  is  strongly  curved)  and,  thus,  one  would  anticipate  the  inlet  conditions  to 
be  of  major  importance.  Recall,  however,  that  all  the  calculations  for  the  1980  Stanford 
Conference  started  with  the  measured  inlet  flow  conditions.  As  a  result  (see  figure  28a)  almost  all 
the  predictions  were  quite  satisfactory  at  0  =  0®  but  major  disagreements  were  observed  at  0  =  90”. 
Interestingly  enough,  at  the  exit  of  the  bend,  most  of  the  Stanford  Conference  calculations  exhibit 
the  same  general  features  as  the  present  results.  For  instance,  most  of  the  methods  predicted 
reasonably  well  the  radial  velocity  component,  but  did  not  predict  as  high  streamwise  mean 
velocities  in  the  core  as  indicated  by  the  experimental  data  (only  one  method  predicted  velocities  of 
the  correct  magnitude  but  then  the  predicted  high  speed  core  was  displaced  close  to  the  outer  wall). 
The  common  factor  between  the  present  numerical  method  and  most  of  the  methods  used  in  the 
1980  Stanford  Conference  is  the  use  of  the  isotropic  k-e  model.  We  should  point  out,  of  course, 
that  in  the  present  method  the  calculations  are  carried  out  all  the  way  to  the  wall,  while  all  the 
Stanford  Conference  methods  employed  the  wall-function  approach  to  bridge  the  gap  between  the 
logarithmic  layer  and  the  wall.  For  the  flow  under  consideration,  high  levels  of  turbulence 
ani.sotropy  are  present  throughout  the  bend.  Near  the  inlet  of  the  bend,  for  instance,  the  stress- 
driven  secondary  motion,  generated  in  the  straight  entry  tangent,  affects  the  dynamics  of  the  flow 
development  but  only  locally,  since  it  is  been  quicklv  taken  over  by  the  pressure-driven  secondary 
motion.  On  the  other  hand,  near  the  exit  of  the  bend,  Humphrey  et  al.  (1983)  reported  high  levels 
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of  turbulence  anisotropy,  caused  by  the  redistribution  of  the  Reynolds  stresses  by  the  secondary 
motion,  as  well  as  regions  of  negative  turbulent  kinetic  energy  production  near  the  inner  convex 
wall.  These  phenomena  cannot  be  accounted  for  by  an  isotropic  eddy  viscosity  model  and, 
therefore,  the  failure  to  accurately  simulate  this  flow  must  be  attributed  largely  to  the  turbulence 
model. 

IX.6  180"  square  bend  of  Chang  (1983) 

This  test  case  corresponds  to  the  measurements  of  Chang  (1983)  who  carried  out  detailed 
mean  and  turbulent  flow  measurements  through  a  strongly-curved  180"  square  bend  and  its 
downstream  tangent  at  Re  =  56,700  and  De  =  21,900.  An  inlet  tangent  of  31  hydraulic  diameters 
long  was  used  in  order  to  ensure  nearly  fully-developed  flow  at  the  entrance  to  the  bend. 

The  calculations  are  carried  out  on  a  numerical  grid  with  119x41x22  nodes  in  the 
streamwise,  radial  and  normal  directions,  respectively.  The  streamwise  spacing  inside  the  bend  is 
2"  while  the  first  coordinate  surface  just  off  the  duct  walls  is  located  at  y"^  =  2  almost  everywhere. 
The  computational  domain  starts  5  hydraulic  diameters  upstream  of  the  inlet  to  the  bend  and 
extends  7  diameters  downstream  of  the  end  of  the  bend.  The  inlet  distributions  of  the  mean 
velocity  and  turbulence  quantities  are  specified  from  a  straight  duct  calculation  as  in  the  previous 
test  case. 

The  calculated  velocity  profiles  at  several  streamwise  locations  are  compared  with  the 
measurements  in  figure  32.  The  discrepancies  between  experiment  and  calculations  observed  at  the 
upstream  station  -Idh  are  due  to  the  inaccuracies  in  the  inlet  conditions  caused  by  the  inability  of 
the  turbulence  model  to  predict  the  stress-driven  secondary  motion.  At  the  next  downstream 
station,  9  =  3",  the  calculations  are  in  good  overall  agreement  with  the  measurements,  although 
discrepancies  associated  with  the  prediction  of  the  boundary  layer  thickness  on  the  concave  wall, 
the  reasons  for  which  have  already  been  discussed,  can  be  clearly  seen,  particularly  near  the 
bottom  wall  of  the  duct  (see  profiles  at  D  and  E).  Further  downstream,  however,  large 
discrepancies  between  experiment  and  computations,  similar  to  those  encountered  in  previous 
turbulent  flow  calculations,  appear  in  the  region  between  the  duct  centerline  and  the  inner  wall. 
More  specifically,  at  9  =  90"  the  pressure-driven  secondary  motion  transports  fluid  from  the 
boundary  layer  of  the  convex  wall  towards  the  duct  centerplane,  resulting  in  the  big  "holes"  in  the 
streamwise  velocity  observed  in  the  A,  B  and  C  profiles.  The  calculations  fail  completely  to 
reproduce  this  feature  of  the  flow  field,  a  trend  which,  once  more,  indicates  that  the  strength  of  the 
secondary  motion  is  grossly  underpredicted.  Similar  discrepancies  between  experiment  and 
calculations,  although  not  as  severe  due,  perhaps,  to  the  reduction  of  the  strength  of  the  secondary 
motion,  persist  at  the  next  station,  9  =  130". 
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In  a  recent  study,  Choi  et  al.  (1989)  have  calculated  the  same  bend  geometry  using  a 
anisotropic,  algebraic  Reynolds-stress  closure  and  a  computational  grid  of  similar  size  as  the  one 
used  in  the  present  calculations.  Their  calculations  were  certainly  more  successful  than  the  present 
ones  in  predicting  at  least  qualitatively  the  shapes  of  the  streamwise  velocity  profiles,  although 
significant  quantitative  differences  remained  between  the  experiment  and  their  computation.  They 
also  reported  a  multi-cellular  structure  of  the  calculated  secondary  motion  at  the  9  =  130^^  station 
with  the  main  vortex  breaking  down  into  three  smaller  vortices.  No  definite  conclusion  can  be 
drawn,  however,  about  how  successful  their  calculations  were  as  far  as  the  prediction  of  the 
secondary  motion  is  concerned  since  no  detailed  comparisons  with  the  data  were  reported. 

X.  FLOW  IN  STRAIGHT  AND  CURVED  TRANSITION  DUCTS 

In  this  section  we  report  a  series  of  laminar  and  turbulent  flow  calculations  with  typical 
transition  duct  configurations.  Two  test  cases  are  studied:  i)  turbulent  flow  through  a  straight 
circular-to-rectangular  transition  duct  for  which  detailed  turbulent  flow  measurements  have  been 
recently  reported  by  Davis  and  Gessner  ( 1992),  and  ii)  laminar  and  turbulent  flow  through  a  typical 
hydroturbine  draft  tube.  The  draft  tube  is  a  strongly  curved  diffuser  whose  cross-sectional  shape 
also  changes  from  circular  at  the  inlet  to  rectangular  at  the  exit.  No  measurements  are  available  for 
the  latter  geometry.  Although  some  measurements  in  a  similar  geometry  have  been  reported,  it  has 
not  been  possible  to  acquire  the  necessary  geometrical  details  to  carry  out  a  meaningful  calculation. 

X.l  Turbulent  flow  in  a  straight  circular-to-rectangular  transition  duct 

The  measurements  of  Davis  and  Gessner  (1992)  were  carried  out  at  a  Reynolds  number. 
Re  =  3.9x10^,  based  on  the  bulk  velocity  and  the  inlet  diameter.  The  circular-to-rectangular  (CR) 
duct  configuration  chosen  for  their  experiment  has  an  exit  aspect  ratio  of  3.0  and  a  transition  length 
of  1.5  inlet  diameters  over  which  changes  in  the  cross-sectional  shape  occur.  At  each  streamwise 
location,  the  cross-sectional  shape  is  defined  by  the  equation  of  a  super-ellipse.  Details  about  the 
precise  geometry  definition  can  be  found  in  Davis  (1992).  The  geometry  of  the  CR  duct  along 
with  the  stations  where  measurements  were  carried  out  are  shown  in  figure  33. 

The  computational  domain  extends  from  station  1  to  station  6  and  a  grid  topology  similar  to 
that  used  for  pipe  bend  calculations  is  employed.  The  numerical  grid,  generated  algebraically  using 
the  EAGLE  grid  generation  code  (Thompson,  1987),  consists  of  46x51x27  grid  nodes  in  the  axial, 
radial  and  tangential  directions,  respectively.  The  grid  lines  are  concentrated  near  the  duct  wall 
using  hyperbolic  tangent  stretching  functions.  The  first  coordinate  surface  just  off  the  solid  wall  is 
located,  almost  everywhere,  at  y+  =  0.5  with  8  to  10  points  within  the  sublayer  and  the  buffer 
layer.  Typical  cross-sectional  views  of  the  numerical  mesh,  as  well  as  a  view  of  the  mesh  on  the 
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duct  wall,  are  given  in  figure  34.  Due  to  the  symmetry  of  the  duct  and  symmetric  inlet-flow 
conditions,  only  one  quadrant  of  the  duct  is  simulated. 

Inlet  conditions  are  specified  at  station  1  using  the  very  detailed  data  of  Davis  and  Gessner 
(1992);  they  reported  mean  velocity  measurements  up  to  the  edge  of  the  sublayer  and  turbulence 
measurements  up  to  y+  =  200.  Detailed  discussion  of  the  flow  conditions  at  station  1  can  be  found 
in  Davis  (1992).  Here  it  suffices  to  say  that  the  inflow  conditions  correspond  to  developing  pipe 
flow  with  a  boundary  layer  thickness  of  approximately  14  percent  of  the  inlet  diameter  and  a 
friction  velocity  of  0.0406. 

Computed  and  measured  distributions  of  wall  static  pressure  coefficient  around  the 
perimeter  at  stations  3  to  6  are  shown  in  figure  35.  In  this  figure,  S  is  distance  along  the  wall 
along  the  perimeter,  and  Sref  is  one  quarter  of  the  duct  perimeter  at  each  station.  The  calculations 
are  in  fair  agreement  with  the  experiment  since  they  reproduce  all  the  measured  trends.  However, 
the  calculated  pressures  are  somewhat  higher  than  the  measured  ones  at  stations  3  and  4  with  this 
trend  reversing  at  the  exit  of  the  duct  (station  6).  This  indicates  that  the  pressure  gradients  are  not 
being  accurately  resolved. 

The  calculated  contours  of  mean  streamwise  velocity  are  compared  with  the  measurements 
at  stations  3,  4,  5  and  6  in  figure  36.  At  station  3,  both  the  calculated  and  measured  contours 
follow  the  general  shape  of  the  wall  and  are  in  good  agreement  with  each  other.  At  the  next 
downstream  station  (station  4),  the  thickening  of  the  measured  boundary  layer  near  the  shorter  wall 
indicates  that  a  pressure-driven  secondary  motion  (induced  by  the  rapid  geometrical  changes)  starts 
to  develop.  The  calculations,  although  in  fair  agreement  with  the  measurements,  predict  a 
significantly  fuller  boundary  layer  in  the  vicinity  of  the  side  wall.  The  distortion  of  the  measured 
isovels  at  stations  5  and  6  is  due  to  a  secondary  flow  pattern  which  develops  into  a  pair  of  vortices 
along  the  shorter  sidewall.  The  calculations  fail  to  reproduce  the  measured  isovel  shapes,  a  trend 
which,  as  in  previously  reported  calculations  with  ducts  of  regular  cross-section,  indicates  that  the 
predicted  secondary  motion  is  much  weaker  than  the  measured  one. 

The  extent  to  which  the  present  calculations  resolve  the  near-wall  region  is  demonstrated  in 
figure  37,  where  measured  and  computed  velocity  profiles  are  plotted  in  wall  coordinates  at 
stations  5  and  6  along  the  vertical  (Y1 -profiles)  and  horizontal  (Y2-profiles)  planes  of  symmetry  at 
each  station.  The  solid  line  in  these  figures  corresponds  to  the  sublayer  and  logarithmic  portions  of 
the  usual  law  of  the  wall.  It  is  seen  that,  for  all  profiles  shown  in  figure  37,  the  calculations  exhibit 
the  correct  asymptotic  behavior  in  the  sublayer.  Moreover,  the  calculated  profiles  are  in  good 
agreement  with  the  measurements  along  the  vertical  plane  of  symmetry  (Y1 -profiles)  where  the 
boundary  layer  is  thin.  Along  the  horizontal  plane  of  symmetry,  however,  the  calculations  agree 
with  the  experiment  only  in  the  inner  layer,  up  to  about  y*"  =  1000,  hut  significant  discrepancies 
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appear  farther  away  from  the  wall.  These  di-screpancies  underscore  the  failure  of  the  calculations  to 
reproduce  the  region  of  low  streamwise  momentum  within  the  developing  vortical  flow. 

X.2  Laminar  and  turbulent  flow  in  a  typical  hydroturbine  draft  tube 

The  draft  tube  is  a  strongly  curved  diffuser  located  right  beneath  the  turbine  in  a 
hydropower  plant.  Its  role  is  to  deliver  the  exhaust  water  from  the  turbine  to  the  tailwater  basin  at  a 
reduced  speed  in  order  to  recover  part  of  the  velocity  head  that  is  not  recovered  by  the  turbine.  The 
velocity  head  recovered  within  the  draft  tube  represents  a  significant  portion  of  the  total  effective 
head  of  the  turbine  and,  therefore,  the  design  of  the  tube  is  of  crucial  importance  for  the  overall 
efficiency  of  the  hydropower  plant.  A  typical  draft  tube  consists  of  a  short  conical  diffuser 
followed  by  a  strongly  curved  elbow  of  varying  cros.s-.section  and  then  a  rectangular  diffuser 
section.  Its  cross-sectional  shape  changes  continuously  from  circular,  at  the  inlet,  to  elliptical 
within  the  elbow,  and  finally  to  rectangular  at  the  exit.  The  flow  that  enters  the  draft  tube-the 
wake  of  the  turbine  blades— is  turbulent  and  three  dimensional  with  high  levels  of  swirl.  The 
already  complex  inlet  flow  undergoes  additional  straining  as  it  passes  through  the  tube,  induced  by 
the  rapid  area  changes  and  the  very  strong  longitudinal  curvature,  resulting  in  an  extremely 
complicated  shear  flow  with  very  strong  vortical  motions,  which  are  often  accompanied  by  regions 
of  streamwise  flow  reversal  at  off-design  operating  conditions. 

Here  we  report  laminar  and  turbulent  flow  calculations  for  a  typical  draft  tube 
configuration.  Our  objective,  herein,  is  twofold:  i)  to  demonstrate  the  feasibility  of  carrying  out 
turbulent  flow  calculations  all  the  way  to  the  wall  for  complex,  three-dimensional  geometries  of 
practical  interest,  and  ii)  to  identify  areas  upon  which  future  research  efforts  should  concentrate  in 
order  to  improve  the  numerical  and  physical  modelling  of  such  flows.  It  is  important  to  point  out 
that  recent  attempts  to  numerically  simulate  the  flow  in  a  draft  tube  (Vu  and  Shyy  (1990), 
Camarero  et  al.  (1991))  have  been  restricted  to  the  use  of  very  coarse  grids  in  conjunction  with  the 
wall-function  approach  for  the  treatment  of  the  near-wall  flow.  As  a  result,  and  given  the  inherent 
complexity  of  the  flow,  the  so  obtained  solutions  have  provided  only  limited  information  about  the 
structure  of  the  flowfield.  The  calculations  to  be  subsequently  reported,  however,  constitute  the 
first  attempt  to  model  the  details  of  the  flow  through  a  draft  tube  using  fine  numerical  meshes. 

X.2.1  Draft  tube  geometry  and  computational  grid 

The  draft  tube  configuration,  used  for  the  present  computations,  is  based  on  one  of  the 
draft  tubes  at  the  Norris  Power  Plant  in  Tennessee,  which  is  designed  to  operate  with  66.()()0  H.P. 
turbines.  The  geometry  of  the  draft  tube  was  made  available  by  the  Tennessee  Valley  Authority 
(Waldrop,  1991b).  The  neat  lines  of  the  tube  along  with  its  plan  and  elevation  views  are  shown  in 
figures  38  and  39.  The  area  expansion  ratio  for  this  draft  tube  (ratio  of  the  exit  over  the  inlet  cross- 
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sectional  area)  is  approximately  4.5:1  while  the  radius  of  curvature  of  the  elbow  is  1.34  diameters 
of  the  inlet  circular  cross-section.  The  configuration  shown  in  figures  38  and  39  has  two  piers, 
symmetrically  placed  about  the  centerline,  supporting  the  downstream  rectangular  diffuser.  For 
simplicity,  particularly  in  grid  generation,  the  piers  are  omitted  in  the  present  calculations  but  their 
effect  is  accounted  for  by  appropriately  .scaling  the  dimensions  of  the  cross-sections  such  that  the 
longitudinal  variation  of  the  net  flow  area  remains  the  same.  The  various  cro.ss-.sections  of  the 
scaled  draft  tube,  without  the  piers,  are  shown  in  figure  40. 

A  grid  topology  similar  to  that  used  for  the  pipe-bend  calculations  (see  figure  20)  is  also 
adopted  for  the  draft  tube.  The  physical  and  transformed  domains  along  with  the  associated 
coordinate  systems  are  shown  in  figure  41.  Note  that,  for  the  present  applications— unlike  all  our 
previous  pipe  calculations— the  inlet  velocity  profile  is  not  in  general  symmetric  and  thus  the  entire 
tube  geometry  has  to  be  considered.  The  transformation  of  the  physical  to  the  computational 
domain  is  achieved  by  introducing  an  artificial  cut  along  the  symmetry  plane  of  the  tube  and 
mapping  it  into  two  planes  (ABCD  and  HGFE  in  figure  41)  where  periodic  conditions  are  applied. 
With  reference  to  figure  41,  the  overall  computational  domain  consists  of  the  inlet  (ABHG)  and 
exit  (DCFE)  planes,  the  solid  wall  boundary  (BCFG),  the  two  periodic  boundaries  (ABCD  and 
HGFE)  and  the  tube  centerline  singular  boundary  (ADEH). 

The  computational  grid  for  each  cross-section  is  generated  algebraically  using  the  EAGLE 
grid  generation  code  (Thompson,  1987).  The  grid  lines  are  concentrated  near  the  walls  using  the 
hyperbolic  tangent  stretching  function.  As  in  all  previous  calculation^  particular  care  is  exercised 
to  keep  the  maximum  stretching  ratio  near  1.3  everywhere.  The  cross-sectional  grids  are  then 
stacked  along  the  centerline  of  the  tube  to  complete  the  three-dimensional  grid.  A  view  of  the 
three-dimensional  grid  on  the  surface  of  the  draft  tube  is  shown  in  figure  42,  while  typical  cross- 
sectional  views  are  shown  in  figure  43. 

X.2.2  Boundary  conditions 

With  reference  figure  41,  the  boundary  conditions  for  the  calculations  are  applied  as 
follows: 

Inlet:  The  di.stributions  of  the  three  velocity  components,  and  the  turbulence  parameters  (k, 
e)  in  the  case  of  turbulent  flow,  are  specified  at  the  inlet  plane  (ABGH).  For  the  laminar  flow 
calculations,  a  plug  flow  profile  with  zero  swirl  is  specified  for  the  velocity  field.  For  the  turbulent 
calculations,  several  velocity  profiles  with  and  without  inlet  swirl  are  considered  (see  subsequent 
sections).  The  pressure  is  computed  from  the  interior  nodes  using  linear  extrapolation. 

Periodic  boundaries:  The  governing  equations  are  solved  on  the  periodic  boundaries 
(planes  ABCD  and  HGFE)  in  a  similar  fashion  as  for  any  interior  node  by  appropriately 
introducing  fictitious  periodic  lines. 
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Solid  wall:  No-slip,  no-flux  conditions  for  the  velocity  components,  k  =  0,  and  linear 
extrapolation  for  the  pressure  are  employed. 

Exit  boundary:  To  facilitate  the  application  of  outflow  boundary  conditions  the 
computational  domain  is  extended  by  adding  a  straight  tangent  (3  to  4  inlet  diameters  long)  at  the 
end  of  the  draft  tube.  At  the  end  of  this  straight  tangent,  the  velocity  components  and  the 
turbulence  quantities  are  computed  by  assuming  zero  streamwise  diffusion  while  the  pressure  is 
computed  using  linear  extrapolation.  The  downstream  extension  of  the  computational  domain  was 
found  necessary  in  order  to  ensure  that  no  reverse  flow  reaches  the  exit  boundary,  in  which  case 
the  calculations  become  unstable  and  eventually  fail  to  converge. 

X.2.3  Laminar  flow 

The  laminar  flow  calculations  are  carried  out  at  Reynolds  number  Re  =  1,000,  based  on  the 
bulk  velocity  and  the  diameter  of  the  circular  inlet  cross-section.  A  plug  flow  profile  without  swirl 
is  specified  at  the  inlet.  Since  the  inlet  profile  is  symmetric  only  half  of  the  draft  tube  needs  to  be 
computed.  Two  numerical  grids  are  used  for  the  calculations  in  order  to  study  the  grid  dependency 
of  the  computed  solutions:  i)  a  coarse  grid  with  43x47x51  (13  planes  in  the  elbow)  nodes,  and  ii) 
and  a  fine  grid  with  64x47x51  (31  planes  in  the  elbow)  nodes  in  the  streamwise,  tangential  and 
radial  directions,  respectively.  The  minimum  grid  spacing,  just  off  the  solid  wall,  is  4xl0'-^  inlet 
diameters  for  both  grids.  Note  that  the  near-wall  resolution  of  the  two  grids  is  the  same  as  that 
used  for  the  subsequently  reported  turbulent  flow  calculations  in  order  to  test  the  robustness  of  the 
numerical  methoo  on  highly  stretched  grids. 

The  solutions  obtained  with  the  coarse  and  fine  grids  are  compared  in  figures  44  and  45. 
In  figure  44,  the  streamwise  velocity  profiles  at  the  end  of  the  elbow  (section  XV  in  figure  38)  are 
plotted  along  the  Horizontal  and  vertical  cross-sectional  axes  of  symmetry,  while  figure  45  depicts 
the  contours  of  constant  static  pressure  on  the  symmetry  plane.  The  agreement  between  the  coarse 
and  fine  grid  solutions  is  good,  since  only  relatively  small  discrepancies  can  be  ob.served.  All  the 
subsequently  reported  results  have  been  obtained  on  the  fine  grid. 

The  comp  ited  velocity  vectors  on  the  symmetry  plane  of  the  draft  tube,  along  with  a  blow 
up  of  the  downstream  region,  are  shown  in  figure  46.  The  flow  in  the  first  half  of  the  elbow 
exhibits  all  the  typical  features  of  a  plug  flow  entering  a  strongly  curved  passage.  More 
specifically,  the  flow  is  accelerated  near  the  inner  (top)  side  of  the  tube  since  it  is  driven  by  the, 
initially,  favorable  longitudinal  pressure  gradient  Along  the  outer  side  (bottom  wall),  on  the  other 
hand,  the  initially  adverse  longitudinal  pressure  gradient  induces  a  small  region  of  reversed  flow. 
Further  downstream,  however,  a  strong  secondary  motion  develops— driven  by  the  transverse 
pressure  gradients— which  transports  high  momentum  fluid  towards  the  outer  (bottom)  wall  of  the 
elbow  (see  figure  48).  As  a  result,  the  flow  is  accelerated  along  the  outer  wall  while  it  is  retarded 
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along  the  inner  (top)  wall.  This  process,  in  conjunction  with  the  continuous  increase  of  the  cross- 
sectional  area  in  the  streamwise  direction,  leads  to  a  large  region  of  reversed  flow  at  the  top  of  the 
draft  tube.  This  originates  near  the  end  of  the  elbow  and  extends  all  the  way  to  the  exit  from  the 
tube. 

Figures  47  and  48  show  the  contours  of  streamwise  velocity  and  the  cross-flow  vectors  at 
sections  IX,  XV  and  XXXI  (see  figure  41).  At  section  IX  (located  at  9  =  45®)  a  very  strong 
secondary  motion  develops  with  two  distinct  swirls  near  the  left  and  right  side  walls.  The 
secondary  motion  tends  to  sweep  the  flow  away  from  the  vertical  axis  of  symmetry  and  as  a  result 
the  maxima  of  the  streamwise  velocity  appear  in  the  vicinity  of  the  side  walls.  At  the  exit  of  the 
elbow  (section  XV)  the  magnitude  of  the  secondary  motion  is  significantly  reduced  and  the  two 
swirls  shift  towards  the  corner  between  the  inner  (top)  and  end  walls.  A  region  of  reversed  flow 
appears  in  the  vicinity  of  the  top  wall  which  induces  the  acceleration  of  the  flow  observed  near  the 
outer  wall.  At  the  exit  of  the  draft  tube  (section  XXXI)  the  secondary  motion— although  very 
weak— sweeps  the  high  speed  fluid  towards  the  two  side  walls  while  a  large  separated  region  is  still 
present  near  the  top  wall. 

X.2.4  Turbulent  flow 

A  series  of  calculations  was  carried  out  at  a  typical  model-scale  Reynolds  number.  Re  = 
1 .  IxlO^  (based  on  the  bulk  velocity  and  the  inlet  diameter  of  the  draft  tube),  in  order  to  investigate 
the  effect  of  the  inlet  velocity  profile  on  the  flow  development  through  the  draft  tube.  Two 
different  inlet  conditions  were  simulated:  i)  a  fully-developed  pipe-flow  profile  (obtained  from  a 
straight  pipe  computation),  combined  with  a  free-vortex  with  an  axial-velocity  defect  (case  FV), 
which  approximates  the  real  flow  situation  since  it  accounts  for  the  wake- like  outflow  from  the 
turbine  runner;  and  ii)  the  fully-developed  pipe-flow  profile  for  the  streamwise  velocity  combined 
with  a  solid-body  rotation  swirl .  In  both  cases,  the  swirl  velocity  at  the  wall  is  reduced  to  zero 
using  the  inner  leg  of  Johnston's  (1960)  triangular  model  relating  the  secondary  and  primary 
velocity  components  in  a  three-dimensional  turbulent  boundary  layer.  The  swirl  intensity  S  is 
defined,  in  the  present  study,  as  the  ratio  of  the  area  average  swirl  velocity  at  the  inlet  divided  by 
the  bulk  velocity  (for  .solid  body  rotation  this  definition  yields  S  =  0.66Vm,  where  V^  is  the 
maximum  swirl  velocity).  According  to  this  definition,  solid-body  swirl  intensities  of  S  =  0,  0.33 
and  0.66  (cases  S1,S2,S3,  respectively)  were  investigated  in  the  present  study.  The  inlet  velocity 
profiles  for  all  the  cases  computed  are  shown  in  figure  49.  In  the  following  sections  we  discuss 
the  grid  sensitivity  of  the  computed  solutions  and  compare  the  computed  flow  structures 
corresponding  to  different  inlet  velocity  profiles  and  swirl  intensities. 
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Grid  dependence  study 

The  grid  dependence  study  was  performed  only  for  the  case  FV.  since  this  corresponds  to 
the  real  life  flow  situation  through  the  draft  tube.  The  calculations  were  carried  out  on  two 
numerical  grids,  a  coarse  one  with  45x94x5 1  nodes,  and  a  fine  one  with  64x94x5 1  nodes  in  the 
streamwise,  circumferential  and  radial  directions,  respectively.  The  streamwise  spacing  inside  the 
elbow  is  7.5®  for  grid  A  and  3.75®  for  grid  B.  For  both  grids,  the  first  coordinate  surface  just  off 
the  solid  wall  is  located  everywhere  such  that  1  <  y+  <  4,  with  at  least  4  coordinate  surfaces  placed 
within  the  sublayer  and  the  buffer  layer.  Approximately  4000  iterations  are  required  for  three 
orders  of  magnitude  reduction  of  the  velocity  and  pressure  residuals,  which,  for  the  fine  grid,  take 
about  4.5  hours  of  CRAY-YMP  CPU  time. 

The  solutions  computed  on  the  coarse  and  fine  grids  are  compared  in  figures  50  to  53.  The 
piezometric  pressure  coefficient  on  the  wall,  around  the  perimeter  at  the  end  of  the  elbow  (section 
XV),  is  plotted  in  figure  5^  i  i  both  grids.  Here,  S  is  the  distance  along  the  wall  and  Sref  is  the 
total  perimeter  at  that  station.  It  is  seen  that  both  solutions  exhibit  similar  trends,  but  the  finer  grid 
predicts  higher  overall  pressures  as  well  as  a  steeper  pressure  gradient  along  the  inner  wall  of  the 
bend  (the  midpoint  of  the  inner  wall  is  located  at  S/Sref  =  0.5).  The  contours  of  piezometric 
pressure  in  the  plane  symmetry  of  the  draft  tube  are  shown  in  figure  5 1 .  The  general  features  of 

both  solutions  are  very  similar  in  the  region  between  the  inlet  and  the  beginning  of  the  elbow. 

* 

Note,  for  instance,  the  sharp  rise  of  pressure  in  the  inlet  conical  diffuser  as  indicated  by  the  closely 
spaced  pressure  contours—in  fact  most  of  the  pressure  recovery  appears  to  take  place  in  this  inlet 
region.  Further  downstream,  however,  and  towards  the  end  of  the  elbow,  the  fine  grid  solution 
predicts  steeper  streamwise  pressure  gradients  as  well  as  higher  pressures  near  the  outer  wall. 

In  figure  52  the  streamwise  velocity  profiles,  computed  on  both  grids,  are  plotted  along  the 
horizontal  and  vertical  axes  of  symmetry  at  the  end  of  the  elbow.  Along  the  horizontal  line  of 
symmetry  (figure  52a)  both  solutions  are  in  close  agreement,  although  the  fine  grid  solution 
predicts  somewhat  lower  velocities  near  the  left  endwall.  Significant  discrepancies  between  the 
two  computed  solutions  appear,  however,  along  the  vertical  axis  of  symmetry  (figure  52b)  and 
particularly  in  the  vicinity  of  the  inner  wall  where  the  coarse  grid  solution  predicts  a  much  thicker 
boundary  layer  as  compared  to  the  fine  grid  prediction.  A  more  global  picture  of  the  computed 
solutions  is  given  in  figure  53,  where  the  contours  of  constant  streamwise  velocity  are  plotted  at 
the  end  of  the  elbow  and  the  end  of  the  draft  tube  (sections  XV  and  XXXI).  The  general  structure 
of  the  isovels  appears  to  be  the  similar  for  both  grids  at  both  sections  but  a  clo.ser  look  reveals 
significant  differences.  At  the  end  of  the  elbow  (figure  53a),  for  instance,  the  closely  spaced 
isovels  near  the  inner  wall,  computed  on  the  fine  grid,  indicate  fuller  velocity  profiles  in  that  region 
(compare,  for  example,  the  0.5  and  0.6  isovels  in  both  solutions).  This  trend  implies  that  the  fine 
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grid  solution  predicts  a  stronger  secondary  motion  than  the  corresponding  prediction  on  the  coarse 
grid,  and  this  tends  to  sweep  the  isovels  closer  to  the  inner  wall,  resulting  in  the  fuller  velocity 
profiles  in  that  region.  At  the  exit  of  the  draft  tube  (figure  53b),  the  coarse  grid  solution  predicts 
lower  velocities  in  the  center  of  the  cross-section  with  regions  of  reversed  flow  both  on  the  inner 
and  outer  walls.  The  fine  grid  solution,  on  the  other  hand,  predicts  only  one  small  region  of 
reversed  flow  near  the  inner  wall. 

Finally,  figure  54  depicts  the  calculated  velocity  vectors  on  the  plane  of  symmetry  as  well 
as  a  blow  up  of  the  downstream  diffuser  region.  The  two  solutions  exhibit  similar  overall  features 
through  the  elbow  and  in  the  downstream  diffuser.  The  coarse  grid  solution,  however,  predicts  a 
much  larger  region  of  reversed  flow  near  the  inner  wall  of  the  downstream  diffuser  which  is 
accompanied— for  continuity  to  be  satisfied— by  higher  velocities  in  the  vicinity  of  the  outer  wall. 
This  trend  is  consistent  with  the  discrepancies  between  the  two  solutions  already  observed  in  the 
streamwise  velocity  contours  and  profiles.  Also,  the  reduction  of  the  separated  flow  region  in  the 
fine  grid  solution  is  consistent  with  the  steeper  streamwise  pressure  gradients,  observed  in  the 
pressure  contours  computed  on  the  fine  grid  (figure  51). 

The  above  comparisons  clearly  indicate  that,  unlike  the  laminar  flow  calculations  where  the 
coarse  and  fine  grid  solutions  were  quite  similar,  the  streamwise  resolution  of  the  coarse  grid  is  not 
sufficient  to  capture  the  evolution  of  the  streamwise  and  secondary  tlow  through  the  elbow  and  in 
the  downstream  diffuser  (recall  that  the  coarse  grid  has  only  13  planes  within  the  elbow).  The 
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coarse  grid  solution,  however,  exhibits  all  the  general  features  observed  on  the  fine  grid  and  can, 
thus,  be  used  to  predict  quantitatively  the  flow  field.  For  that  reason  and  in  order  to  save  computer 
resources,  all  the  subsequently  reported  calculations  with  different  inflow  conditions  are  carried  out 
on  the  coarse  grid.  Note  that  the  main  objective  of  these  calculations  is  to  get  an  overall  idea  about 
the  impact  that  the  inflow  conditions  have  on  the  flow  development  through  the  draft  tube  rather 
than  to  perform  detailed  comparisons  with  experimental  data,  which  are  not  available  at  this  time. 

Effect  of  inflow  conditions  on  the  flow  development 

The  velocity  distributions  in  the  symmetry  plane  of  the  draft  tube  computed  with  the  four 
different  inflow  conditions  are  shown  in  figure  55.  The  figure  also  includes  an  enlargement  of  the 
exit  region  of  the  draft  tube  so  that  the  effects  of  inflow  conditions  can  be  easily  identified.  For 
case  S 1  (no  inlet  swirl,  figure  55a),  the  maximum  of  the  velocity  profile  initially  shifts  towards  the 
inner  wall  of  the  elbow  (due  to  the  initially  favorable  pressure  gradient  along  this  wall)  but  this 
trend  is  reversed  half  way  through  the  elbow  by  the  strong  pressure-driven  secondary  motion. 
Note  that  this  behavior  is  consistent  with  experimental  data  for  curved  ducts  of  regular  cross- 
section.  Towards  the  exit  from  the  elbow,  the  strong  adverse  pressure  gradient,  induced  by  the 
continuous  area  increase,  results  in  a  large  recirculation  or  reverse  axial-flow  zone  covering  the 
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area  between  the  centerline  and  the  inner  wall,  with  reverse-flow  velocities  as  high  as  30  percent  of 
the  bulk  velocity.  The  favorable  effect  of  a  solid-body  swirl  on  the  flow  development  can  be  seen 
in  the  subsequent  figures  55b  and  55c.  For  case  S2,  for  instance,  the  extent  of  the  recirculation 
zone  in  the  downstream  diffuser  is  significantly  reduced  and  the  maximum  backflow  velocities  do 
not  exceed  15  percent  of  the  bulk  velocity,  while  for  case  S3  the  recirculation  zone  disappears 
entirely.  For  this  latter  case  it  is  interesting  to  note  that,  starting  half  way  through  the  elbow,  the 
flow  bifurcates  towards  the  inner  and  outer  walls  and,  as  a  result,  the  minimum  of  the  velocity 
appears  near  the  centerline  of  the  draft  tube.  For  the  FV  case,  shown  in  figure  55d,  a  region  of 
very  low  velocity  appears  near  the  centerline  (due  to  a  combination  of  the  adverse  pressure  gradient 
and  the  wake-like  inlet  velocity  profile)  and  the  maximum  of  the  velocity  remains  near  the  inner 
wall  throughout  the  elbow.  In  the  downstream  diffuser,  a  small  recirculation  zone  also  appears 
near  the  inner  wall  but  the  velocity  profiles  between  the  centerline  and  the  outer  wall  tend  to  be 
more  uniform  as  compared  to  the  corresponding  profiles  of  cases  S 1  and  S2. 

The  effect  of  the  inflow  conditions  on  the  pressure  field  is  depicted  in  figure  56,  where 
contours  of  constant  piezometric  pressure  are  plotted  in  the  symmetry  plane  for  all  four  cases.  The 
pressure  in  these  figures  has  been  referenced  with  respect  to  the  pressure  at  the  center  of  the  inlet 
cross-section.  For  cases  SI,  S2  and  S3  it  is  seen  that  the  overall  pressure  rise  between  the  inlet  of 
the  draft  tube  and  the  end  of  the  elbow  increases  with  increasing  swirl  intensities.  In  addition,  for 
case  S3,  almost  50  percent  of  the  overall  pressure  rise  takes  place  within  the  inlet  conical  diffuser. 
For  the  FV  case  the  pressure  rises  very  rapidly  in  the  inlet  region  and  almost  70  percent  of  the 
overall  pressure  rise  is  achieved  within  the  inlet  conical  diffuser. 

The  contours  of  streamwise  velocity  and  the  associated  secondary  flow  vectors  at  stations 
IX,  XV  and  XXXI  are  shown  in  figures  57  to  59  (the  inner  wall  of  the  draft  tube  is  at  the  top).  It 
can  be  seen  that  the  inlet  swirl  has  a  dramatic  effect  on  the  flowfield  through  the  elbow  and  in  the 
downstream  diffuser.  At  station  IX,  the  two  curvature-induced  swirls  located  near  the  inner  (top) 
wall,  symmetrically  about  the  vertical  centerplane  in  case  SI  (figure  57a),  transport  high 
momentum  fluid  towards  the  outer  (bottom)  wall  along  the  centerplane,  and  this  results  in  the  low 
velocity  regions  observed  near  the  inner  wall.  The  imposition  of  an  inlet  swirl  which  counteracts 
the  pressure-driven  secondary  motion,  however,  tends  to  reverse  this  trend.  In  case  S2.  for 
instance  (figure  57b),  the  maximum  of  the  streamwise  velocity  has  shifted  near  the  inner  wall 
while,  for  the  S3  case  (figure  57c),  the  maximum  velocity  also  appears  near  the  inner  wall  but  a 
region  of  high  velocity  has  appeared  near  the  outer  wall  as  well.  This  accumulation  of  high 
momentum  fluid  near  the  inner  wall  increases  the  resistance  of  the  flow  against  the  adverse 
pressure  gradient  and,  thus,  results  in  the  reduction  of  the  recirculation  zone  in  the  downstream 
diffuser  as  discussed  above.  A  similar  trend  is  also  obvServed  in  case  FV  (figure  57d). 
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At  the  end  of  the  elbow  (figure  58),  a  large  region  of  reversed  flow  appears  near  the  inner 
wall  in  case  the  SI.  As  the  swirl  intensity  increases,  however,  the  separated  flow  region  is 
significantly  reduced  and  eventually  eliminated  for  sufficient  high  swirl  intensities  (cases  S3  and 
FV).  For  case  S2,  the  maximum  of  the  velocity  is  located  in  the  vicinity  of  the  inner  wall  to  the 
right  of  the  vertical  axis  of  symmetry.  For  the  S3  case,  on  the  other  hand,  the  bifurcation  of  the 
high  speed  flow  towards  the  inner  and  outer  walls,  observed  in  figure  55c,  persists  in  this  station 
as  well.  In  addition,  a  core  of  low  speed  fluid  is  seen  to  develop  near  the  junction  between  the 
inner  wall  and  the  left  endwall.  For  the  FV  case,  the  maximum  of  the  velocity  is  also  located  near 
the  inner  wall  and  to  the  right  but  now  a  larger  cross-sectional  area  near  the  right  endwall  is 
occupied  by  high  speed  fluid. 

At  the  exit  from  the  draft  tube  (figure  59),  zero  inlet  swirl  (figure  59a)  results  in  a  large 
region  of  reversed  flow,  near  the  inner  wall,  and  in  a  secondary  motion  that  sweeps  the  high  speed 
flow  towards  the  two  end  walls.  For  the  S2  case,  the  separated  flow  region  is  significantly 
reduced  and  a  single  swirl  appears  near  the  right  endwall.  High  velocity  fluid  is  accumulated  near 
both  the  left  and  the  right  endwalls  but  the  maximum  velocity  occurs  in  the  vicinity  of  the  right 
endwall.  Further  increase  of  the  swirl  intensity  (case  S3)  results  in  an  accumulation  of  high  speed 
fluid  almost  exclusively  near  the  right  endwall  with  a  .small  region  of  reversed  flow  at  the  junction 
between  the  inner  and  left  end  walls.  A  flow  pattern  somewliat  similar  to  that  of  case  S2  is  also 
observed  at  the  exit  of  the  tube  in  case  FV  (figure  59d).  More  specifically,  high  speed  fluid  is 
found  in  the  vicinity  of  both  endwalls  with  the  maximum  velocity  near  the  right  endwall. 
However,  a  region  of  low  velocity,  with  two  separation  bubbles  near  the  inner  and  outer  walls,  is 
found  in  the  middle  of  the  cross- section. 

Finally,  the  near-wall  resolution  achieved  by  the  present  numerical  method  and  turbulence 
model  is  demonstrated  in  figure  60,  where  velocity  profiles  are  plotted  in  wall  coordinates  for  cases 
SI,  S2  and  S3,  and  compared  with  the  universal  law  of  the  wall.  The  profiles  shown  in  these 
figures  are  plotted  at  two  axial  locations,  at  the  exit  from  the  elbow  and  the  at  the  end  of  the  draft 
tube,  on  the  symmetry  plane  from  the  outer  wall  towards  the  centerline.  At  the  exit  from  the 
elbow,  all  three  velocity  profiles  collapse  into  a  single  curve  and  conform  with  the  law  of  the  wall 
up  to  y*"  <  80,  approximately.  This  behavior  is  consistent  with  the  conclusions  of  the  recent 
experimental  work  of  Davis  and  Gessner  (1992).  Their  data  indicate  that  at  the  end  of  the 
transition  duct  the  measured  velocity  profiles  exhibit  logarithmic  behavior  for  values  of  y*"  less  than 
80,  but  they  fall  below  the  logarithmic  curve  for  greater  values  of  y*".  At  the  end  of  the  draft  tube, 
however,  the  computed  velocity  profiles  conform  with  the  law  of  the  wall  only  within  the  sublayer 
but  none  of  them  appears  to  exhibit  a  logarithmic  region.  More  specifically,  the  velocity  profile  for 
S 1  (no  swirl)  falls  well  below  the  logarithmic  curve  while  increased  swirl  intensities  tend  to  rever.se 
that  trend.  This  latter  behavior  is  consistent  with  the  .strong  adverse  pressure  gradient  which  the 
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flow  experiences  in  the  downstream  diffuser.  The  behavior  of  the  near- wall  flow  for  case  SI,  on 
the  other  hand,  can  be  attributed  to  the  strong  divergence  of  the  flow  away  from  the  centerplane 
towards  the  left  and  right  endwalls  (a  divergence  induced  by  the  pressure-driven  secondary 
motion).  As  the  inlet  swirl  intensity  increases  it  counteracts  the  secondary  motion  and,  as  a  result, 
the  effects  of  the  adverse  pressure  gradient  dominate  in  the  near-wall  region.  The  velocity  plots 
clearly  show  the  location  of  the  near-wall  grid  points  and  the  resolution  that  is  achieved  in  the 
sublayer. 


XI.  SUMMARY  AND  CONCLUSIONS 

Two  numerical  methods  were  described  and  their  results  compared  with  each  other  in  order 
to  identify  the  relative  merits  and  disadvantages  of  these  methods  when  applied  to  calculate  laminar 
and  turbulent  flows  in  curved  ducts  with  varying  cross-section.  Method  1  employs  the  finite- 
analytic  discretization  for  the  convective  and  viscous  terms  and  solves  the  fully-transformed 
governing  equations  in  generalized  curvilinear  coordinates  using  the  ADI  method.  Method  II,  on 
the  other  hand,  is  a  finite-difference  method  which  integrates  in  time  the  partially-transformed 
governing  equations  (with  the  cartesian  velocity  components  retained  as  unknowns  in  generalized 
curvilinear  coordinates)  using  the  four-stage,  explicit  Runge-Kutta  algorithm,  enhanced  with  local 
time  stepping  and  implicit  residual  smoothing.  For  turbulent  flows,  both  methods  employ  the  two- 
layer,  two  equation  k-e  model  of  Chen  and  Patel  (1988).  The  methods  were  applied— on  the  same 
numerical  grid,  with  identical  boundary  conditions,  starting  from  the  same  initial  conditions,  and 
using  the  same  convergence  criteria— to  calculate  laminar  and  turbulent  flows  through  strongly 
curved  ducts  of  regular  cross-section.  The  computed  solutions  were  compared  with  each  other  and 
with  experimental  data  in  order  to  assess  the  spatial  resolution  of  each  method. 

Method  II  was  subsequently  applied  to  calculate  laminar  and  turbulent  flows  in  several 
curved  ducts  of  square  and  circular  cross-sections,  as  well  as  turbulent  flow  in  a  straight  circular- 
to-rectangular  transition  duct,  and  a  curved  ciicular-to-rectangular  transition  duct,  namely  a 
hydroturbine  draft  tube.  The  objective  of  these  calculations  was  twofold:  i)  to  further  validate  the 
spatial  resolution  and  numerical  performance  of  Method  II;  and  ii)  to  validate  the  ability  of  the  two- 
layer  turbulence  model  to  predict  the  origin,  growth  and  decay  of  pressure-driven  vortical 
structures  in  complex,  three-dimensional  shear  flows.  To  facilitate  the  validation  process,  a 
detailed  literature  survey  was  carried  out  in  order  to  summarize  ail  the  experimental  and 
computational  work  in  the  area  of  duct  flows  and  identify,  among  the  available  experimental 
studies,  those  that  are  more  appropriate  for  validating  numerical  methods  and  turbulence  closures. 
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The  study  is  concluded  with  application  of  Method  II  to  simulate  the  flow  through  a  typical 
hydroturbine  draft  tube.  Unlike  previous  attempts  to  simulate  such  flows,  the  present  method 
resolves  the  near- wall  flow  using  fine  computational  meshes  across  the  sublayer.  Numerical 
solutions  were  obtained,  at  typical  model-scale  Reynolds  numbers  for  various  inlet  swirl  intensities 
and  axial  velocity  profiles,  for  a  draft-tube  configuration  based  on  the  installation  at  the  Norris  Dam 
in  Tennessee. 

The  following  general  conclusions  are  drawn  from  the  present  study: 

i)  An  accurate  simulation  of  complex,  three-dimensional  flows  requires  very  careful 
discretization  of  the  governing  flow  equations— particularly  the  convective  terms— in  order  to  ensure 
that  numerical  diffusion  is  kept  to  a  minimum  and  real  viscous-flow  processes  are  not  obscured. 
This  was  demonstrated  in  the  comparison  of  the  two  numerical  methods,  since  the  finite-analytic 
method  (Method  I)  appeared  to  consistently  underpredict  the  strength  of  the  secondary  motion  and 
consequently  its  impact  on  the  streamwise  flow  development.  This  is  attributed  to  the  fact  that  the 
finite-analytic  discretization  automatically  switches  to  a  first-order  accurate  upwind  di.scretization  in 
regions  of  the  flow  where  the  cell  Reynolds  number  exceeds  a  certain  limiting  value.  The 
performance  of  Method  I  in  highly  three-dimensional  flows  considered  here  could  be  improved  if 
the  hybrid  nine-point  discretization  formula  employed  in  the  present  version  of  the  method  is 
replaced  with  a  more  accurate  discretization  involving  additional  points  in  the  upstream  and 
downstream  sections  of  the  finite-analytic  element.  This,  however,  would  result  in  additional 
computational  effort.  Method  n,  which  uses  a  .second-order  accurate  upwind  scheme,  was  shown 
to  be  sufficiently  accurate  to  predict-at  least  for  the  laminar  flow  cases  where  no  turbulence-model 
related  uncertainties  are  present— the  evolution  of  the  streamwise  and  the  secondary  flow  with 
remarkable  accuracy. 

ii)  The  present  study  indicates  that  isotropic  eddy-viscosity  based  turbulence  models  are 
inadequate  for  resolving  the  origin,  growth  and  decay  of  pressure-driven  secondary  motion  in 
complex,  three-dimensional  shear  flows,  even  when  the  most  advanced,  state-of-the-art  near-wall 
models  are  used  to  resolve  the  flow  all  the  way  to  the  wall.  In  all  the  turbulent  flow  calculations 
performed  herein,  the  two-layer  k-e  model  failed  to  predict  the  strength  of  the  secondary  motion  as 
well  as  its  effect  on  the  streamwise  flow  development.  These  di.screpancies  cannot  be  attributed  to 
numerical  discretization  errors  since  the  calculations  were  carried  out  using  Method  II,  the 
numerical  performance  of  which  was  carefully  established  through  a  series  of  complex  laminar- 
flow  calculations.  While  it  is  well  known  that  turbulence  models  of  the  type  employed  here  cannot 
predict  Reynolds-stress-driven  secondary  motion,  the  reasons  for  the  failure  of  the  two-layer 
model  to  capture  the  pressure-driven  secondary  motion  are  not  entirely  clear,  particularly  because  it 
appears  to  resolve  the  important  near-wall  layers  in  even  the  most  complex  flow  situations.  One 
may  only  speculate  that  the  observed  discrepancies  are  due,  at  least  in  part,  to  the  equally  well 
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known  deficiency  of  the  k-e  model  to  account  for  the  increased  levels  of  turbulence  production 
along  concave  walls.  Whether  it  is  this  factor,  or  the  general  anisotropy  of  the  Reynolds  stresses, 
cannot  be  answered  without  further  computations  using  turbulence  models  that  claim  greater 
generality.  This  is  obviously  a  topic  for  further  research.  Extensions  of  Method  II  to  incorporate 
models  based  on  the  individual  Reynolds-stress  transport  equations  are  now  under  way. 

iii)  The  computations  with  the  draft  tube  geometry  demonstrated,  for  the  first  time,  the 
feasibility  of  carrying  out  turbulent  flow  calculations  all  the  way  to  the  wall  for  complex,  three- 
dimensional  configurations.  The  numerical  method  was  shown  to  yield  robust  and  efficient 
solutions  on  highly  stretched,  very  fine  computational  grids— with  as  many  as  310.(){)()  grid  nodes. 
The  efficiency  of  the  numerical  method  can  be  significantly  improved,  however,  by  implementing 
the  multigrid  acceleration  technique  which  can  reduce  the  overall  CPU  time  required  for 
convergence  by  a  factor  of  four  to  five. 

iv)  The  inlet  swirl  intensity  and  the  axial  velocity  profile  were  found  to  have  a  dramatic 
effect  on  the  flow  development  throughout  the  draft  tube.  The  numerical  results  indicate  that  high 
swirl  intensities  tend  to  diminish  the  stalling  characteristics  of  the  flow  through  the  elbow  and  to 
increase  the  overall  pressure  recovery.  In  addition,  the  comp  ..ff  olutions  were  shown  to 
reproduce  with  remarkable  accuracy  the  correct  limiting  behavior  of  the  flow  within  the  sublayer. 
While  the  present  solutions  are  qualitatively  similar  to  those  reported  by  others  for  different  draft- 
tube  geometries,  no  specific  conclusions  can  bv  irawn  at  this  time  insofar  as  the  overall  accuracy 
of  the  computed  solutions  is  concerned  since  no  experimental  data  is  available  to  validate  the 
numerical  predictions.  An  experiment  is  under  way  at  Voith  Hydro  Inc.  to  measure  the  flow  in  a 
scaled  model  of  the  Norris  Dam  draft  lube,  including  the  flow  dividing  piers  in  the  downstream 
reach  of  the  tube.  The  numerical  method  will  be  modified  to  handle  such  a  geometry,  and 
comparisons  will  be  made  between  the  calculations  and  data  as  they  are  gathered. 
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Table  1.  Experiments  in  laminar  flow 
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length  of  downstream  tangent 
bend's  turning  angle  in  degrees 


Table  2.  Experiments  in  turbulent  flow  in  curved  ducts 
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Table  3.  Experiments  in  turbulent  flow  in  straight  transition  ducts 
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AR  :  aspect  ratio  at  the  exit  station 
C-R  :  circular-to-rectangular  transition  duct 
R-C  :  rectangular-to-circular  transition  duct 
S-C  :  square-to-circular  transition  duct 
Ltr  :  length  of  transition 


Typical  flow  structure  in  the  coma*  region  (from  Gessner  and  Jones,  1965) 
(a)  isotach  pattern;  (b)  secondary  flow  pattern 


Fig.  3  Hnite  difference  computational  ceil 


91 


Fig.  4  Curvilinear  coordinates  and  computational  grid  for  a  curved  square  duct 

(a)  physical  solution  domain;  (b)  transformed  solution  domain 
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Fig.  4  Cundlinear  coordinates  and  computational  grid  for  a  curved  square  duct 
(c)  Typical  cross-sectional  and  plane  of  symmetry  grid  views 


93 


IduaL ) 


\ 


Fig.  5  Comparison  of  the  convergence  histories  for  Methods  I  and  11 

(Dtict  of  Humphrey  et  al.  (1977);  Re=790;  Grid  A) 
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Q.O 


Fig.  7  \feasuTed  and  calculated  streamwise  velocity  profiles  along  radial  lines; 

o,  measurements  (Humphrey  ct  ai.,  1977); - ,  Method  I;  ,  Method  n 

(a)  z  =  •0.25dh;  (b)  z  =  0.0 
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Measured  and  calculated  streamwise  velcKity  profiles  along  z-lines: 
o,  measurements  (Humphrey  et  al.,  1977); - ,  Method  I; - ,  Method  II 
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Fig.  9  Calculated  radial  velocity  profiles  along  z-lines  (duct  of  Humphrey  et  al.,  1977): 
- ,  Method  I; - ,  Method  II 
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Fig,  10  Measured  (Humphrey  et  al,,  1977)  and  calculated  streamwise  velocity  contours 
(a)  0  =  600 
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Calculation 

Method  I 

Experiment 


Calculation 

Method  H 
(b)  0  =  90° 

Fig.  10  Measured  (Humphrey  et  al.,  1977)  and  calculated  streamwise  velocity  contours 
(b)  0=900 
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Experiment 


Calculation 

Method  I 
Experiment 


Calculation 

Method  n 
(a)  Xj^=-0.25  d^ 

Fig.  1 1  Measured  (Taylor  et  al.,  1982)  and  calculated  streamwise  velocity  contours 
(a)  xh  =  -  0.25dh 
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Experiment 


Calculation 

Method  I 

Experiment 


Calculation 


Method  n 
(b)  0  =  30® 

Fig.  1 1  Measured  (Taylor  et  al.,  1982)  and  calculated  streamwise  velocity  contours 
(b)0  =  3Oo 
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Experiment 


Calculation 
Method  I 

Experiment 


Calculation 

Method  n 
(c)  0  =  60° 


Fig.  1 1  Measured  (Taylor  et  al.,  1982)  and  calculated  streamwise  velocity  contours 

(c)  0  =  600 
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Experiment 


Calculation 

Method  I 


Method  n 
(d)  0=775'* 


Fig.  1 1  Measured  (Taylor  et  al.,  1982)  and  calculated  streamwise  velocity  contours 
(d)  e  =  77.5° 
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Experiment 


Calculation 

Method  H 
(e)  Xj^=+0.25d^ 


Fig.  1 1  Measured  (Taylor  et  al.,  1982)  and  calculated  streamwise  velocity  contours 
(e)  Xh  *  +  0.2^ 
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Calculation 

Method  I 


Experiment 


Fig.  1 1  Measured  (Taylor  ct  al.,  1982)  and  calculated  streamwise  velocity  conttturs 
(f)  Xh  =  +  2.5^ 


men 
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Measured  and  calculated  radial  velocity  profiles  along  z-lines: 
o,  measurements  (Taylor  et  al.,  1982); - ,  Method  I; - ,  Method  IF 


Top  View 


Fig.  14  Curved- wall  wind  tunnel,  coordinates  and  locations  of  measurement  stations  for 
the  duct  of  Kim  (1991) 
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Fig.  IS  Measured  (Kim,  1991)  top-wall  boundary  layer  at  station  U1 


0.0  0.5  1.0  0.0  0.5 


Fig.  16  Measured  and  calculated  mean  velocity  profiles  along  radial  lines 

□ ,  measurements  (Kim,  1991); - ,  Method  I; - ,  Method  II 

(b)  Streamwise  mean  velocity;  Station  45 
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Fig.  16  Measured  and  calculated  mean  velocity  profiles  along  radial  lines 

□  ,  measurements  (Kim,  1991); - ,  Method  I; - ,  Method  11 

(c)  Streamwise  mean  velocity;  Station  75 


00 


Fig.  16  Measured  and  calculated  mean  velocity  profiles  along  radial  lines 

O,  measurements  (Kim,  1991); - ,  Method  I; - ,  Method  11 

(d)  Streamwise  mean  velocity;  Station  D1 


114 


Fig.  16  Measured  and  calculated  mean  velcKity  profiles  along  radial  lines 

□,  measurements  (Kim,  1991); - ,  Method  I; - ,  Method  II 

(e)  Streamwise  mean  velocity;  Station  D2 
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Fig.  16  Measured  and  calculated  mean  velocity  profiles  along  radial  lines 

□,  measurements  (Kim,  1991); - ,  Method  I; - ,  Method  II 

(i)  Radial  mean  velocity.  Station  D1 


119 


so- 


station  02 
Z-0.25 


Station  02 
Z-0.50 


0.0  0.5  1.0  0.0  0.5  1.0 


Y/H  Y/H 

Fig.  16  Measured  and  calculated  mean  velocity  profiles  along  radial  lines 

□,  measurements  (Kim,  1991); - ,  Method  I; - ,  Method  11 

(j)  Radial  mean  velocity;  Station  D2 


Fig.  16  Measured  and  calculated  mean  velocity  profiles  along  radial  lines 

□,  measurements  (Kim,  1991); - ,  Method  I; - ,  Method  II 

(k)  Vertical  mean  velocity;  Stadon  15 
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Fig.  16  Measured  and  calculated  mean  velocity  profiles  along  radial  lines 

□,  measurements  (Kim,  1991); - ,  Method  I; - ,  Method  II 

(1)  Vertical  mean  velocity;  Station  45 
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Fig.  16  Measured  and  calculated  mean  velocity  profiles  along  radial  lines 

□,  measivements  (Kim,  1991); - ,  Method  I; - ,  Method  II 

(m)  Vertical  mean  velocity.  Station  75 
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0.0  0.5  1.0  0.0  0.5  1.0 
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Fig.  16  Measured  and  calculated  mean  velocity  profiles  along  radial  lines 

□,  measurements  (Kim,  1991); - ,  Method  I; - ,  Method  11 

(n)  Vertical  mean  velocity;  Station  D1 
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Fig.  16  Measured  and  calculated  mean  velocity  pnrfiles  along  radial  lines 

□,  measurements  (Kim,  1991); - ,  Method  I; - ,  Method  11 

(o)  Vertical  mean  velocity;  Station  D2 
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0.0  0.5  1.0  0.0  0.5  1.0 


Y/H  Y/H 

Fig.  17  Measured  and  calculated  turbulent  kinetic  energy  profiles  along  radial  lines 

O,  measurements  (Kim,  1991); - ,  Method  I; - ,  Method  II 

(a)  Station  45 
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Fig.  17  Measured  and  calculated  turbulent  kinedc  energy  profiles  along  radial  lines 

□,  measurements  (Kim,  1991); - ,  Method  I; - ,  Method  II 

(b)  Station  D1 
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station  U2  Station  15  Station  45  Station  75  Station  D1  Station  D2 

Fig.  1 8  Stream  wise  velocity  contours  (U  =1 .09  to  0.99  by  0.02, 0.95  to  0.65  by  0.05) 

(a)  Measurements  (Kim,  1991) 


Stream  wise  velocity  contours  (U  =1.09  to  0.99  by  0.02,  0.95  to  0.65  by  0.05) 
(b)  Method  I 
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Fig.  18  Streamwise  velcx^ity  contours  (U  =1 .09  to  0.99  by  0.02,  0.95  to  0.65  by  0.05) 
(c)  Method  II 


Fig.  20  Curvilinear  coordinates  for  a  typical  pipe  bend  configuradon 
(a)  physical  soludon  domain;  (b)  transformed  soludon  domain 


(a)  d--0° 


Fig.  22  Measured  (Bovendeerd  et  al.,  1987)  and  calculated  streamwise  velocity  conK>urs 
(a)  0=00;  (b)  0=  11.7o 
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Calculation 


(c)  0  =  23.4® 


Fig.  22  Measured  (Bovendeerd  ct  al.,  1987)  and  calculated  streamwise  velocity  contours 
(c)  0  *  23.40;  (d)  0  »  39.8° 


(e)  a=58.5“ 


(f)  0=81.9® 

Fig.  22  Measured  (Bovendeerd  et  al.,  1987)  and  calculated  streamwise  velocity  contours 
(e)  0  =  58.50;  (f)  9  =  81.9° 
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0=81.9° 


0=58.5° 


0=39.8° 


0  =  234° 


0=117° 


0  =  46° 


0=0° 


Fig.  23  Measured  and  calculated  streamwise  velocity  profiles  on  the  plane  of  symmetiy: 
o,  measurements  (Bovendeerd  et  al.,  1987); - .calculation 


Experiment 


(a)  0  =  0° 


Fig.  25  Measured  (Rowe,  1970)  and  calculated  total  velocity  head  contours 
(a)0  =  OO;(b)  0  =  30° 
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Fig.  25  Measured  (Rowe,  1970)  and  calculated  total  velocity  head  contours 
(c)0  =  45O;(d)  0  =  60° 
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Experiment 


(f)  0=120° 


Fig.  25  Measured  (Rowe,  1970)  and  calculaied  total  velocity  head  contours 
(e)  0  =  900;  (0  0  =  120o 
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;  (b)  circumferential  mean  velocity 


Fig.  27  Calculated  mean  secondary  flow  vectors  and  mean  streamwise  velocity  contours 
(pipe  bend  of  Azzola  et  at.,  1984) 

(a)  e  =  3°;  (b)  0  =  45° 


ig.  27  Calculated  mean  secondary  flow  vectors  and  mean  streamwise  velocity  contours 
(pipe  bend  of  Azzola  et  al.,  1984) 


Experiment 


Fig.  29  Calculated  streamwise  velocity  contours,  for  fully  developed  straight  dua  flow, 
compared  with  the  measurements  of  Humphrey  et  al.  (1981)  at  Xh  =  -2.5dh 
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5 


Calculation 
(c)  0=71° 

Experiment 


Calculation 

(d)  0=90° 


Fig.  30  Measured  (Humphrey  et  al.,  1981)  and  calculated  streamwise  velocity  contours 
(c)  0  =  710;  (d)  0  =  900 
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Experiment 


Fig.  3 1  Measured  (Humphrey  et  al.,  1981)  and  calculated  radial  velocity  contours 
at  0  =  90° 
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Measured  and  calculated  streamwise  velocity  profiles  along  radial  lines: 
o,  measurements  (Chang,  1983); - .calculation 


Fig.  32  Measured  and  calculated  strcamwise  velocity  profiles  along  radial  lines: 

o,  measurements  (Chang,  1983); - .calculation 

(c)  e  =  90°;  (d)  e  =  130° 


Fig.  33  Coordinates  and  locations  of  measurements  stations  for  the  transition  duct  of 
Davis  and  Gessner  (1992) 


56 
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Fig.  34  Typical  views  of  the  computational  grid  (CR  duct  of  Davis  and  Gessner,  1992) 
(a)  3-D  surface  view 
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Fig.  34  Typical  views  of  the  computational  grid  (CR  duct  of  Davis  and  Gcssncr,  1992) 
(b)  Cross-sectional  views 
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Fig.  35  Measured  (Davis  and  Gcssncr,  1992)  and  computed  peripheral  wall  pressure 
coefficient  distributions 
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z/R 


Fig.  36  Measured  (Davis  and  Gessner,  1992)  and  computed  axial  velocity  contours 
(a)  Station  3 


2/R 


z/R  2/R 


-0.6  — - - — ^ ^ ^ ^ ^ ^ ! 

0.0  0.5  1.0  1.5 


y/R 

Fig.  36  Measured  (Davis  and  Gessner,  1992)  and  computed  axial  velocity  contours 
(c)  Station  5;  (d)  Station  6 
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Fig.  37  Measured  (Davis  and  Gessner,  1992)  and  computed  law-of-the-wall  velocity 
profiles  along  the  horizontal  and  vertical  planes  of  symmetry 
(a)  Station  5 
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Fig.  37  Measured  (Davis  and  Gessner,  1992)  and  computed  law-of- the- wall  velocity 
profiles  along  the  horizontal  and  vertical  planes  of  symmetry 
(b)  Station  6 


Fig.  38  Neat  lines  of  the  TVA  Norris  Power  Plant  draft  tube 


Boundary 

(b) 

Fig.  41  Curvilinear  coordinates  for  a  typical  draft  tube  configuration 
(a)  physical  solution  domain;  (b)  transformed  sdution  domain 
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Fig.  42  Typical  3-D  surface  view  of  the  computational  grid  for  the  draft  tube  geometry 
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U/Ub 


(a) 


(b) 


Fig,  44  Grid-dependence  study  (Re=  1000):  Stieamwise  velocity  profiles  at  station  XV 
(a)  vertical  plane  of  symmetry;  (b)  horizontal  plane  of  symmetry 
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(b) 


Fig.  46  Velocity  field  on  the  plane  of  symmetry  (Re=1000) 
(a)  Overall  view;  (b)  Exit  region 
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(c) 


Fig.  47  Streamwise  velocity  contours  (Rc=1000) 

(a)  Station  IX;  (b)  Station  XV;  (c)  Station  XXXI 
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Fig.  48  Secondary  flow  vectors  (Rc=l000) 

(a)  Station  IX;  (b)  Station  XV;  (c)  Station  XXXI 


velocity  profiles  at  the  inlet  of  the  draft  tube 
;  ^1;  (c)  Case  S2:  (d)  Case  S3 


(a) 


Fig.  52  Grid-dqwndence  study  (Rc=l.lxlO®;  Case  FV):  Streamwise  velocity  profiles  at 
station  XW.  (a)  Vertic^  plane  of  symmetry,  (b)  Horizontal  plane  of  symmetry 
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Fig.  53  Grid-dependence  smdy  (Re=l.lxlO^,  Case  FV):  Streamwise  velocity  contours 
(a)  Station  XV 


(b) 


Fine  Grid 


Fig.  53  Grid-dependence  study  (Re=  1.1x10*^,  Case  FV):  Streamwise  velocity  contours 
(b)  Station  XXXI 
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(a 


Fig.  54  Grid-dependence  study  (Re=l.lxlO^;  Case  FV);  (a)  Velocity  field  on  the  plane 
<rf  symmetry  (overall  view) 
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Fine  Grid 


54  Grid-rfependence  study  (Re=l,lxlO^;  Case  FV):  (b)  Velocity  field  on  the  plane 
of  symmetry  (exit  region) 
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Fig.  55  Velocity  field  on  the  plane  of  syirjnetiy  for  various  inflow  conditions 
(a)  Case  SI ;  (b)  Case  , S2 
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Fig.  55  Velocity  field  on  the  plane  of  symmetry  for  various  inflow  conditions 
(c)  Case  S3;  (d)  Case  FV 


0.- 


(b)  Cose  S2 


Fig.  56  Static  pressure  contours  on  the  plane  of  symmetry  for  various  inflow  conditions 
(a)  CaU  SI;  (b)  Case  S2 
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Fig,  56  Static  pressure  contours  on  the  plane  of  symmetry  for  various  inflow  conditions 
(c)  Case  S3;  (d)  Case  FV 
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(c)  Case  S3 


0.2 


(d)  Cose  FV 


Fig.  57  Streamwise  velocity  contours  and  secondary  flow  vectors  at  section  DC  for 
various  inflow  conditions:  (c)  case  S3;  (d)  case  FV 
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itreamwise  velocity  contours  and  secondary  flow  vectors  at  section  XV  for 
various  inflow  conditions:  (a)  case  SI;  (b)  case  S2 


Fig.  59  Streamwise  velociiy  contours  and  secondary  flow  vectors  at  section  XXXI  for 
various  inflow  conditions:  (a)  case  Si;  (b)  case  S2 


Fig.  59  Streamwise  velocity  contours  and  secondary  flow  vectors  at  section  XXXI  for 
various  inflow  conditions;  (c)  case  S3;  (d)  case  FV 


loq^^(y^) 


Fig.  60  Law-of-the-wall  velocity  profiles  for  various  inflow  conditions 
(a)  Section  XV;  (b)  Section  XXXI 
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